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We review some aspects of the realistic implementation of the dynamical mean-field method. We 
extend the techniques introduced in Ref. Q]to include the calculations of transport coefficients. The 
approach is illustrated on Lai-^Sr^TiOs material undergoing a density driven Mott transition. 



I. INTRODUCTION 

In recent years understanding of the physics of strongly 
correlated materials has undergone tremendous increase. 
This is in part due to the advances in the theoretical 
treatments of correlations, such as the development of 
dynamical mean- field theory (DMFT) 1 . The great al- 
lure of DMFT is the flexibility of the method and its 
adaptability to different systems as well as the simple 
conceptual picture it allows us to form of the dynam- 
ics of the system. The mean-field nature of the method 
and the fact that the solution maps onto an impurity 
model, many of which have been thoroughly studied in 
the past, means that a great body of previous work can 
be brought to bear on the solution of models of corre- 
lated lattice electrons. This is exemplified by the great 
many numerical methods that can be employed to solve 
the DMFT equations. 

DMFT has been very successful in understanding the 
mechanism of the Mott transition in model Hamiltoni- 
ans. We now understand that the various concentration 
induced phase transitions can be viewed as bifurcation 
of a single functional of the Weiss field. The phase di- 
agram of the one-band Hubbard model, demonstrating 
that there is a first order Mott transition at finite tem- 
peratures is fully established 1 . Furthermore Landau like 
analysis demonstrates that all the qualitative features are 
quite generic at high temperatures 2 . However the low- 
temperature ordered phases, and the quantitative aspects 
of the spectra of specific materials clearly require realistic 
treatment. 

This triggered realistic development of DMFT in the 
last decade which has now reached the stage that we 
can start tackling real materials from an almost ab ini- 
tio approach^, something which in the past have been 
exclusively in the domain of density functional theories. 
We are now starting to see the merger of DMFT and 
such ab initio techniques and consequently the opportu- 
nities for doing real electronic structure calculations for 
strongly correlated materials which so far were not within 



the reach of traditional density functional theories. 

Density functional theory (DFT)& is canonical exam- 
ple of ab initio approach, very successful in predicting 
ground state properties of many systems which are less 
correlated, for example the elemental metals and semi- 
conductors. However it fails in more correlated materials. 
It is unable to predict that any systems is a Mott insula- 
tor in the absence of magnetic order. It is also not able 
to describe correctly strongly correlated metallic state. 
As a matter of principle DFT is a theory of ground state 
or thermodynamic properties at finite temperatures. It's 
Kohn-Sham spectra cannot be rigorously identified with 
the excitation spectra of the system. In weakly correlated 
substances the Kohn-Sham spectra is a good approxima- 
tion to start a perturbative treatment of the one-electron 
spectra using the GW method—. However this approach 
breaks down in strongly correlated situations, because it 
is unable to produce Hubbard bands. In orbitally ordered 
situations the LDA+U method^ produces the Hubbard 
bands, however this method fails to produce quasipar- 
ticle bands and hence it is unable to describe strongly 
correlated metals. Furthermore, once long range order is 
lost the LDA+U method reduces to LDA and hence it 
becomes inappropriate even for Mott insulators. 

Dynamical mean-field theory is the simplest theory 
that is able to describe on the same footing total en- 
ergies and the spectra of correlated electrons even when 
it contains both quasiparticle and Hubbard bands. Com- 
bined with LDA, one then has a theory which reduces to 
a successful method (LDA) in the weak correlations limit. 
In the static limit, one can show — that LDA+U can be 
viewed as a static limit of LDA+DMFT used in con- 
junction with the Hartree-Fock approximation. There- 
fore LDA+U is equivalent to LDA+DMFT+further ap- 
proximations which are the only justified in static ordered 
situations. Up to now, the realistic LDA band structure 
was considered with DMFT for purpose of computing 
one-electron (photocmission) spectra and total energies. 

Following Refs. lll'Sl I"! in this paper we extend this ap- 
proach to computation of transport properties. Many 
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transport studies within DMFT applied to model Hamil- 
tonians have been carried out, and the strengths (non- 
perturbative character) and limitations (absence of ver- 
tex corrections) are well understood. However applica- 
tions to real materials require realistic computations of 
current matrix elements. 

There are two ways in which DMFT can be used to 
understand the physics of real materials. The simplest 
approach, outlined in Ref. O]0 is closely tied to the idea 
of model Hamiltonians. This requires i) methodology 
for deriving of the hopping parameters and the interac- 
tion constants ii) a technique for solving the dynamical 
mean-field equations and iii) an algorithm for evaluat- 
ing the transport function which enters in the equations 
of transport coefficients. The second direction is more 
ambitious and focus on an integration of i) and ii) using 
functional formulations^. 

In this paper we review the first approach. The em- 
phasis here is in illustration of different aspects of the 
modeling which affect the final answer. This is necessary 
to obtain a balanced approach towards materials calcula- 
tions. There are now many impurity solvers, they differ 
in their accuracy and computational cost. In the present 
paper we use two impurity solvers the Hirsch-Fye Quan- 
tum Monte Carlo (QMC) method ^ and symmetrized 
finite-U NCA method (SUNCA)±I comparing them in the 
context of simplified models without the additional com- 
plications of real materials. Instead we use the SUNCA 
method as an impurity solver to compute the transport 
properties and new developments using Lai-^Sr^TiOs 
as an example material—. For other reviews of realistic 
implementations of DMFT and electronic structure see 
Ref. [H 

In the next section [H] we shortly review a basic dy- 
namical mean-field theory concepts and their applica- 
tion to realistic structure calculations. As computation 
of transport parameters requires knowledge of the self- 
energy coming from DMFT calculations which is based 
on impurity solvers we present a short review of two im- 
purity solvers used in the paper in section II I II Theory 
of the transport calculations is given in section llVl Test 
system used for transport calculations, which is doped 
LaTiOs ceramics, and DMFT results are described in 
section [3 Results of dc- transport calculations are pre- 
sented in section IVTl And finally we come to conclusion 
in homonymous section IVIII 



II. DYNAMICAL MEAN-FIELD THEORY 

A. Realistic DMFT formalism 

A central concept in electronic structure theory is the 
/-model Hamiltonian. Conceptually, one starts from the 
full many body problem containing all electrons and then 
proceeds to eliminate some high-energy degrees of free- 
dom. The results is a Hamiltonian containing only a few 
bands. The determination of the model Hamiltonian is a 



difficult problem in itself, which has received a significant 
attention 14 i 15 i 16 i 17 i 18 i 19 i 20 . The Kohn-Sham Hamiltonian 
is a good starting point for the kinetic part of the Hamil- 
tonian and can be conveniently expressed in a basis of 
linear muffin-tin orbitals (LMTO's) — , which need not 
be orthogonal (see Appendix IB]) ). as 
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where i,j are atomic site indexes, m is orbital one, and 
a denotes spin. 

It is well known that LDA severely underestimates 
strong electron interactions between localized d- and /- 
electrons because the exchange interaction is taken into 
account only approximately via the functional of electron 
density. To correct this situation, the LDA Hamiltonian 
can be supplemented with a Coulomb interaction term 
between electrons in the localized orbitals (here we will 
call them a heavy set of orbitals). The largest contribu- 
tion comes from the Coulomb repulsion between electrons 
on the same lattice site that we will approximate by the 
interaction matrix U % of the heavy shell (h) of atom i as 
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where m are orbital and a are spin indexes. In the diago- 
nal density approximation electron - electron interaction 
matrix m , m i can be represented using screened 

Coulomb and exchange vertexes as 
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which are expressed via Slater integrals FW , i = 0, 2,4, 6 
in the standard manner >22>. The Slater integrals can be 
linked to the average intra-atomic repulsion U and ex- 
change J obtained from, e.g., LSDA supercell procedures 
via U = F° and J = (F 2 + F 4 )/14. The ratio F 2 /F 4 is to 
a good accuracy a constant ~ 0.625 for d-electrons. Us- 
ing the Coulomb and exchange matrices we can rewrite 
the interaction term as 

Hint ^ ^ ] U aa /Tli a Tli a ' — — ^ ^ ^mm' ^ima^im' — a 

ioLOt' imm'a 
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where index a — (m, a) combines the orbital and spin 
indexes. This equation also provides definition of the 
interaction matrix U aa > which will be used further in the 
paper. 

The LDA Hamiltonian already contains a part of the 
local interaction which has to be subtracted to avoid the 
double counting. The full Hamiltonian is thus approxi- 
mated by 



H — Hlda — Hdc + H int = H + 



(5) 
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where H° is the one-particle part of the Hamiltonian 
and will play a role of the kinetic term within a DMFT 
approach. The double counting correction can not be 
rigorously derived within LDA+DMFT. Instead, it is 
commonly assumed to have a simple static Hartree-Fock 
form, just shifting the energies of the heavy set 



(6) 



Here, r is atomic index in the elementary unit cell. The 
simplest approximation commonly used for E^ c is 
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is the total number of electrons in 



where n h = J2 ma 
the heavy shell (see Appendix [Bj). 

In the spirit of DMFT, the self-energy is assumed to be 
local, i.e. fc-independent, and non-zero only in the block 
of heavy orbitals. Therefore it is convenient to partition 
the Hamiltonian and the Green's function into the light 
and heavy set (denoted by I and h, respectively) as 
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where [•••] means matrix inversion, /j, is the chemical 
potential and O is the overlap matrix (see Appendix lA"|l . 

The main postulate of the Dynamical Mean-Field The- 
ory (DMFT) 1 formalism is that the self-energy is local, 
i.e. it does not depend on momentum, Yl(k,u>) — T,(u>). 
This postulate can be shown to be exact in the limit 
of infinite dimensions provided that the hopping param- 
eters between different sites are scaled appropriately. 
Within this approach, the original lattice problem can 
be mapped onto an Anderson impurity model where the 
local Green's function and the self-energy, G\ oc and E, 
are identified with the corresponding functions for the 
impurity model, i.e. 
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(uj) = E(u>) and G imp (w) = Gi oc (lu). (9) 



Equations © along with the trivial identity 
G loc (u,)=J2 G ( k ^)> 



(10) 



constitute a closed set of self-consistent equations. The 
only thing that remains is to solve the Anderson impurity 
model. 

Notice that statement that the self-energy is diagonal 
is the basis dependent statement and if E(z<x> n ) is momen- 
tum independent in one basis and Uk is a unitary trans- 
formation from one basis to another, and LMTO Hamil- 
tonian, Hlda in the new basis is given by UicHldaUI, 
then the self-energy in the new basis E' = L^E («<;,,,) f/J 
is momentum dependent. Therefore DMFT approxima- 
tion, if at all valid, is valid in one basis^. Hence, we will 



work in a very localized basis where the DMFT approx- 
imation is most justified. 

In DMFT we construct the self-energy, E, as a solution 
of an Anderson impurity model with a non-interacting 
propagator (Weiss function) Go 

Si mP = E 4Wft^(^)«vM (ii) 

OtOi' ,TT f 

aa f £h 

where a and a' are running over indexes ma. The Weiss 
function can be linked to the lattice quantities through 
the local Green's function and self-energy being are re- 
lated to each other by the Dyson equation 

Glo^iUnY 1 = GoiiuJny 1 - E(lW n ). (12) 

Combining Eq. JSJ , (fTUfl and lfP2)l we finally obtain 



G 1 (iu n ) = 



1 



(iuj n + n)O k 



-£(iw„). 
(13) 



One can solve a very general impurity model defined by 
the action (|1 If) and Weiss field 1)13(1 . But it is much 
cheaper to eliminate the light (weakly interacting) bands 
and define an effective action in the subspace of heavy 
bands only. In this way, the local problem can be sub- 
stantially simplified. 



1. Downfolding 

When a group of bands is well separated from the oth- 
ers it is clear that a reduced description of the problem is 
possible. In the one-electron approach it goes under the 
name downfolding—. There are several prescriptions to 
carry out this procedure in the context of DMFT. Per- 
haps, the simplest approach is to reduce the Hamiltonian 
(H Q +H int ) (see Eqs. (0, {Q) to the one having Hubbard 
like form for the bands in question. To estimate the hop- 
ping elements one can perform a tight-binding fit. The 
value of U to be used is then reduced, since ones com- 
puted in the constrained density functional calculations, 
it is screened by the bands which have been eliminated. 
After this procedure we arrive to a Hubbard Hamiltonian 
with a small number of bands 



H= £ t 



where i,j run over lattice sites and m, m! (a, a') label the 
orbital (spin) indices of the heavy set of orbitals. 

We can also perform the downfolding at the level of 
DMFT. The starting point is Eq. © which is used to get 
the heavy block Green's function 
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where the quantity is defined as 



Mt 



{u + p)0 7 (fc) - H°(k), 



(15) 



here 7 is a double index which is combination of h and /. 

The low-energy form of Eq. (|14fl can be found by ex- 
panding around zero frequency and to linear accuracy in 
uj we obtain 



G hh {io) = [Zk 1 ^ - H(k) - Z hh 

k 



(16) 



where renormalization amplitude Zk and effective Hami- 
tonian are given by 
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Further, we can choose a new base in the heavy block 
such that the local Green's function has the usual form. 
First we define the average renormalization amplitude 
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(18) 



and choose the transformation matrix such that 

S^Z- 1 S=1. (19) 

The new overlap matrix and effective Hamiltonian be- 
come 

O eff (k) = tfZ^S, (20) 
H eff (k) = S^H(k)S + f ,O eff (k) , (21) 
£ = S^ZhhS, (22) 

and finally the local Green's function in the new base 
takes the form 



G hh = [(<" + M)O e //(fc) - H eff (k) - £]" 

k 

with the Dyson equation 

Ghh = ^Qhh ~ S i 

and corresponding Weiss function defined as 
Gohh = [(w + /i) - A] -1 , 



(23) 



(24) 



(25) 



where hybridization function A regularly behaves at in- 
finity. 

The self-energy of the reduced model l)23|) is local, 
therefore the DMFT treatment is applicable. In this case, 
only the heavy bands need to be considered in calculation 
which greatly simplifies the complexity of the problem. 

The reduced model, however, has in general more com- 
plicated Coulomb interaction matrix than we chose in 



the original model. If we assume, that the diagonal com- 
ponents of Z are much larger than the off-diagonal, we 
obtain the same simple Hubbard-type interaction term 
within reduced model for the heavy block. The Coulomb 
interaction is however screened by the light bands and is 
reduced to 



U a j3 — ► U a f}Z a Zj3. 



(26) 



We can estimate the magnitude of the reduction of 
Coulomb repulsion by evaluating the above quantity ex- 
plicitly. In the case the original base is orthogonal, i.e. 
O = 1, we get 
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(27) 



where N is the number of heavy bands. 

In the case when the reduced model corresponding to 
Eq. (|23|l is degenerate and the basis is orthogonal (the 
overlap matrix O e f / (k) is equal to unity) the self-energy 
and local Green's function are also degenerate and diag- 
onal. Then the momentum sum in Eq. (12.'}[| can be re- 
placed by the integral over energy and the local Green's 
function can be calculated using the standard Hilbcrt 
transformation 



G(iu> n ) 



ds 



D(e) 



(28) 



Here, the density of states D(e) is the density of states of 
the kinetic part of the reduced model H e ff(k) in Eq. 12311 . 

Other groups have used for D(e) rescaled partial DOS 
as discussed in Appendix^! Possibility to use the Hilbcrt 
transformation is substantially simplifies the calculation 
procedure and brings a number of conceptual simplifica- 
tions — . 



2. Upfolding 

Upfolding is a procedure which is "inverse" to the 
downfolding one. One needs to use Eq. to trans- 
form the self-energy obtained from DMFT calculations, 
S, to the one, S^h, which is inserted to the original 
LDA Hamiltonian in order to compute the local Green's 
function (GF) G(ito n ). The local GF with upfolded self- 
energy reads 



G(iujr, 



dk 



(iu n + (jt) 



Hhh Hhi 
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(29) 



where fi is the LDA chemical potential and Hd c is the 
double counting term Eq. (jSJ. Instead of using formula 
Q) , we rather deduce the constant shift of the heavy set 
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of bands by equating the total number of electrons to the 
integral of the spectral function 

fe a/3 

multiplied by the Fermi function. 



III. IMPURITY SOLVERS 

Among many methods used to solve the impurity prob- 
lem we chose the Quantum Monte Carlo method 10 and 
symmetrized finite-U NCA methodii. In this section we 
briefly describe both of them. 



3. Algorithm to solve DMFT equations 

To close the set of DMFT equations, a method to solve 
the local problem is required. In the following, we will 
focus our attention on two impurity solvers: QMC and 
SUNCA. In section ITTTl we will briefly review both meth- 
ods while a detailed comparison between the results ob- 
tained by those two approaches is given in section IV Dl 
Bellow, we summarize basic steps in the DMFT self- 
consistent scheme that delivers the local self-energy - cru- 
cial quantity to calculate transport and optical properties 
of a solid. 

We started the iteration by a guess for the Weiss field 
Gq 1 from which the local Green's function Gi oc was cal- 
culated by one of the impurity solvers. The self-energy 
was then obtained by the use of the Dyson equation 
(|24[) . Momentum summation over the Brillouin zone, Eq. 
(|23|l . delivers a new guess for the local Green's function 
and through the Dyson equation also for the Weiss field 
Qq ■ The iteration is continued until the convergence is 
found to the desired level. The scheme can be illustrated 
by the following flow-chart 

„_! IMP solver „ DE „ DMFT SCC „_i 

y — > g — > — > y , 

where "DE" stands for the Dyson equation l|24l) and 
"DMFT SCC" means the DMFT self-consistent condi- 
tion Eq. or J2HJ- 

The QMC impurity solver is defined in imaginary time, 
r, therefore the following additional Fourier transforma- 
tions between imaginary time and Matsubara frequency 
points are necessary 



IFT 



Go(t) 



QMC 



FT 



G(t) iAG(iw) 



Here FT and IFT are Fourier and inverse Fourier trans- 
formations, respectively. After the self-consistency is 
reached, the analytic continuation is required to obtain 
the real-frequency self-energy. This issue is addressed in 
section ITlI CI 

The SUNCA method is implemented on real frequency 
axis to avoid the ill-posed problem of analytic continu- 
ation. As an input, it requires the bath spectral func- 
tion A c (u>) = — ^IirQq 1 and delivers the local spec- 
tral function Ad(u>) = — ilmG(w) 



A c {w 



SUNCA 



KK 



A{u>) ^G(w). 



The real part of the local Green's function is obtained by 
the use of the Kramers-Kronig relation (KK) . 



A. The Quantum Monte Carlo method 

There are well known advantages and disadvantages of 
the QMC method and our choice is spurred by the fact 
that despite being slower than other methods the QMC 
is well controlled, exact method. As an input the QMC 
procedure gets Weiss function Gq(t) and as an output 
it produces Green's function G(r). We remind reader 
major steps taken for the QMC procedure. Usually one 
starts with impurity effective action S 

Seff = - [ drdr' J2 <£(t)Qo-Ht,t')c*(t') (30) 
1 f 13 

a, a' 

where {c, c + } operators are fermionic annihilation and 
creation operators of the lattice problem, a — {m, a}. 

The first what we should do with the action (|31|) is to 
discretize it in imaginary time space with time step At 
such that /3 = LAt, and L is the number of time intervals 

Seff - E c + a {T)g - l (ry)c a (T') (31) 

q,tt' 

+ 77 ^ U aa 'n a (T)n a ,(T). 



The next step is to get rid of the interaction term U 
by substituting it by summation over Ising-like auxiliary 
fields. The decoupling procedure is called the Hubbard- 
Stratonovich transformation 25 - 26 



exp{— AT{U aa >n a n a > 



1, 

2 {n ° 



i E exp{A QQ 'S , QQ '(r 



V')}} = 



(32) 



Saa'=±l 



where coshA QQ ' = exp( 



At(/„ 



Saa'iri) are auxiliary 



Ising fields at each time slice. 

In the one-band Anderson impurity model we have 
only one auxiliary Ising field S(ti) — ±1 at each time 
slice, whereas in the multiorbital case number of auxil- 
iary fields is equal to number of a, a' pairs, i.e. a C2- Ap- 
plying the Hubbard-Stratonovich transformation at each 
time slice we bring the action to the quadratic form with 
the partition function 



Z = Tr 



Il dctG M 



(33) 
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where GF in terms of auxiliary fields G a 1 reads as 

G:: |Sm ,,(^') = Q^ l {r,r')e v - (e v 1)5 T , T ,, (34) 
with interaction matrix 

V" — ^ ^aa'S aa >(T)(7 aa i , (35) 

where 

fT oa i — +1 for a < a' 
&au' = — 1 for a > a'. 

Once the quadratic form is obtained one can apply 
Wick's theorem at each time slice and make the Gaus- 
sian integration by Grassmann variables to get the full 
interacting GF 

G a (r,r') = ^Tr {Saal} (36) 
G a , { S aal }(r,r')l[detG-} {Saai} . 

a' 

To evaluate summation Eq. I|3t)|) one uses Monte Carlo 
stochastic sampling. The product of determinants is in- 
terpreted as the stochastic weight and auxiliary spin con- 
figurations are generated by a Markov process with prob- 
ability proportional to their statistical weight. More rig- 
orous derivation can be find elsewhere—^. 

Since the QMC method produces results in complex 
time (G(r m ) with r TO = mAr, m = 1...L) and the DMFT 
self-consistency equations make use of the frequency de- 
pendent Green's functions and self-energies we must have 
an accurate method to compute Fourier transforms from 
the time to frequency domain. This is done by rep- 
resenting the functions in the time domain by a cubic 
splincd functions which should go through original points 
with condition of continuous second derivatives imposed. 
Once we know cubic spline coefficients we can compute 
the Fourier transformation of the splined functions ana- 
lytically (see Appendices ICl and |D"|) . 

B. The SUNCA 

In this section we briefly review the second method, 
used to solve the multiorbital Anderson impurity model, 
called Symmetrized finite-U NCA method 11 . The 
method is based on the self-consistent perturbation the- 
ory with respect to the hybridization strength between 
the effective bath and the local system and is there- 
fore exact in the atomic limit. However, this method 
sums up an infinite class of skeleton diagrams and takes 
into account a subclass of singular vertex corrections 
that are necessary to obtain the correct dynamic low- 
energy Fermi-liquid scaled and correct position of the 
Abrikosov-Suhl resonance. The SUNCA approximation 
does not contain only non-crossing diagrams but rather 
all the three-point vertex corrections of ladder-type and 



should not be confused with the usual finite-U non- 
crossing approximation. 

The SUNCA approach presents the advantage to pro- 
vide directly real frequencies one-particle Green's func- 
tions. Local quantities like densities of states can there- 
fore be computed for any regime of parameters without 
having to perform analytical continuation. Furthermore, 
the SUNCA method can be applied to arbitrary multi- 
band degenerate Anderson impurity model with no ad- 
ditional numerical cost. This is an important advantage 
compared to some other methods like Quantum Monte 
Carlo or exact diagonalization. The method is especially 
relevant for systems with large orbital degeneracy such 
as systems with /-electrons. 

The pathologies that severely limit the usefulness of 
the non-crossing approximation in the context of DMFT 
are greatly reduced with inclusion of ladder-type vertex 
corrections. Nevertheless, they do not completely re- 
move the spurious peak that forms at temperatures sub- 
stantially below the Kondo temperature. To overcome 
this shortcoming, we employed an approximate scheme 
to smoothly continue the solution down to zero tempera- 
ture. This was possible because at the break-down tem- 
perature the solution of SUNCA equations shows an on- 
set of the Fermi-liquid state. As we will show in the 
subsequent chapters by comparison with QMC, SUNCA 
gives correct quasiparticle renormalization amplitude Z 
and the real part of the self-energy at zero frequency ap- 
proaches the Luttinger value. The imaginary part of the 
self-energy, however, has a narrow spurious deep on top 
of the parabola that is formed around zero frequency. To 
remove the deep, we matched the Fermi-liquid parabolic 
form for imaginary part of the self-energy in the small- 
window of the deep such that it smoothly connected the 
intermediate frequency region where the parabola was 
formed. We numerically found that this SUNCA pathol- 
ogy is rapidly reduced with increasing the number of 
bands i.e. it is much less severe in the case of three-band 
model than in one-band case. 

Next, we give some details of the auxiliary particle 
technique together with the definition of the SUNCA 
approximation. Within the auxiliary diagrammatic 
method, the local degrees of freedom can be represented 
by auxiliary particles. To each eigenstate of the local 
Hamiltonian H\ oc |n) = E n \n) we assign an auxiliary op- 
erator a n such that \n) = a^ n \vac). A general Anderson 
impurity model can then be expressed by 

H = E n a\a n + 2J £k m c\ m Ck m + 

n k.m 

^Fn^W'Ckm + h.c), (37) 

where F™ n , — {n\d) m \n'), d\ n is a creation operator for 

the local electron and c\ m creates an electron in the bath 
and m stands for the spin and band index. In order 
that electrons are faithfully represented by the auxiliary 
particles, two conditions must be satisfied: 
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(i) An auxiliary particle a n must be boson (fermion) if 
the state \n) contains even (odd) number of electrons. 

(ii) The local charge Q = J2 n a n a n must be equal to one 
at all times Q — 1 expressing the completeness relation 
for the local states ^2 n \n) (n| = 1. 

The first condition merely request some care that has to 
be taken in evaluating diagrams while the second con- 
straint, projection onto the physical Hilbert space, is 
somewhat more involved but can still be done exactly. 
The term XQ can be added to the Hamiltonian and the 
limit A — * oo has to be taken after the analytic contin- 
uation to the real frequency axes is performed. Taking 
this limit, actually leads to a substantial simplification 
of the analytic continuation. Namely, when evaluating 
self-energies for the pseudo-particles, only the integrals 
around the branch-cuts of the bath electron Green's func- 
tion have to be considered while the integrals around the 
auxiliary particle Green's functions vanish by the projec- 
tion. 

The physical local Green's function (electron Green's 
function in Q=l subspace) can eventually be calculated 
with the help of the Abrikosov trickSl which states that 
the average of any local operator that vanishes in the 
Q = subspace is proportional to the grand-canonical 
(all Q values allowed) average of the same operator 



(38) 



By realizing that the local Green's function is propor- 
tional to the bath electron T-matrix, we have 



G 



1 



loc 



& V 3 (Q>, 



-s CJ 



(39) 



where E c is the bath electron self-energy calculated in 
the grand-canonical ensemble. 

In the case of degenerate Anderson impurity model, 
an important simplification occurs in the Hamiltonian 
Eq. I|37|). Namely, pseudo-particles, corresponding to the 
states with the same number of electrons on the impu- 
rity site, are degenerate with energies Em = —M/i + 
UM(M — l)/2, and pseudo Green's functions 

G M (ioj) = l/(tw - A - E M — SmM), (40) 

where M is the number of electrons on the impurity. In 
the case of Anderson impurity model with N/2 bands, 
we need to consider only N + 1 different propagators in- 
stead of dealing with 2 N pseudo-particles. Furthermore, 
if U is large and we are interested in doping levels not 
too far from an integer filling with M electrons, i.e. close 
to the Mott-insulating state, it is reasonable to assume 
that only fluctuations between local states with M — 1, 
M and M + 1 electrons on the impurity need to be con- 
sidered. The single particle spectra will thus consist of 
lower Hubbard band, upper Hubbard band and a quasi- 
particle resonance while we ignore other Hubbard bands 
that are even more far away from the chemical potential. 




" km r bf 




" km r af 



\f 



FIG. 1: The two bare vertices considered in the case of 
degenerate Anderson impurity model. Throughout this pa- 
per, conduction-electron c propagator is represented by solid 
line, while wiggly, dashed and zigzag lines correspond to the 
pseudo-particles with M — 1, M and M + 1 electrons on the 
impurity, respectively. 



For example, in the case the filling is less than 1.5 we take 
into account empty state, singly occupied and doubly oc- 
cupied states on the impurity while the triple occupancy 
is neglected. This is equivalent to adding to the original 
Hamiltonian an infinite three-particle interaction. 

In the following, we will use letter b to denote the aux- 
iliary particles with M — 1 electrons on the impurity, / 
for M and a for M + 1 electrons on the local level, respec- 
tively. In Fig. ^ we show the two bare three-point ver- 
tices that are left to us in the case of degenerate Ander- 
son impurity model when considering those three types 

( N 

of states. Note that pseudo-particle b is 



de- 



generated, / is 



N 



and a is 



degenerated, 



v M - 1 

N 

M J " i0 \ M+l 
respectively. In case M is odd (even), particles a and b 
are bosons (fermions) while / are fermions (bosons). 

The SUNCA approximation is a conserving approxi- 
mation defined by a Luttinger-Ward type functional $ 
from which all self-energies are obtained as a functional 
derivatives, S a = ^t-. The building blocks of $ are 
dressed Green's functions of pseudo-particles G& (de- 
picted as a wiggly line), G/ (dashed line), G a (zigzag 
line) and bath electron G c (solid line). Due to the exact 
projection, only pseudo-particles are fully dressed while 
bath electron Green's function is not dressed. 

The choice of diagrams was motivated by the 
Schrieffer-Wolf transformation showing that both fluc- 
tuations from M to M — 1 and M to M + 1 have to be 
considered totally symmetrically. In the case of one-band 
model we have for the exchange coupling J = J\ + J-j = 
V 2 (j^ + jj^ji )■ The exact Kondo temperature is propor- 
tional to exp(— tv /J) = exp(— 7r/(Ji + h))- It is well 
known that if one takes into account only non-crossing 
diagrams, the Kondo scale of the resulting approxima- 
tion steeply drops with decreasing U and is of the order 
of exp(— TV I J\ — tv / J2) which can be orders of magnitude 
wrong. This happens because the simultaneous fluctua- 
tions between all three types of pseudo-particles are ne- 
glected. As was shown in Ref. llll one needs to sum up an 




+ 



FIG. 2: Diagrammatic representation of the SUNCA generat- 
ing functional to describe the degenerate Anderson impurity 
model. 



infinite number of skeleton diagrams (SUNCA diagrams) 
to recover the correct exchange coupling consisting of two 
terms J\ and J2 , generated from the two types of vertices 
depicted in Fig. ^ 

The SUNCA Luttinger-Ward functional, shown in 
Fig. consists, in addition to NCA contributions (first 
two diagrams), of the diagrams where two conduction 
electron lines cross only ones and the rest of the electron 
lines cross twice. Diagrams not included in SUNCA have 
higher order crossings. The lowest order neglected are of 
CTMA-typeS&Si where all conduction electrons cross ex- 
actly twice. Note that due to the projection, any contri- 
bution to the Luttinger-Ward functional consists of a sin- 
gle ring of pseudo-particles since at any moment in time 
there must be exactly one pseudo-particle in the system. 

The SUNCA self-energies, obtained by differentiating 
the Luttinger-Ward functional, are shown in Fig. [3J The 
part of the pseudo-particle ring, where conduction elec- 
trons cross exactly twice, can be rewritten in the ladder 
type T-matrix. The additional two conduction lines, that 
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FIG. 3: Diagrammatic representation of the self-energies, 
derived from SUNCA Luttinger-Ward functional (Fig. in 
terms of the renormalized hybridization vertices, defined in 
Fig- El In eacn l me the third diagram is subtracted in order to 
avoid double counting of terms within the first two diagrams. 




FIG. 4: Diagrammatic representation of the Bethe-Salpeter 
equations for ladder-type three-point vertex function. 



cross only once, enable us to close T-matrix and convert 
it into three-point vertex shown in Fig. 0] This is im- 
portant from practical point of view since the tree-point 
vertex is numerically much more tractable than the full 
four-point vertex. 

The expressions for the self-energies £ M defined by 
Fig. and Fig. together with the definitions of the 
Green's functions, Eq. l]4Up. constitute a set of non-linear 
integral equations for = a,b,f. The local Green's 

function is calculated from the Eq. I|39f) using the self- 
consistcntly determined auxiliary propagators G M . The 
SUNCA equations are given explicitly in Appendix [E] 

The SUNCA equations have been evaluated numeri- 
cally by iteration. The two self-consistent loops, SUNCA 
and DMFT, can be merged together into one single set 
of equations. Starting with the initial guess for the 
bath spectral function A c and pseudo-particle Green's 
functions G M , the first guess to the T-matrix defined 
in Fig. 0| is determined. With this, the pseudo-particle 
self-energies as well as the local Green's function may 
be deduced. From the DMFT self-consistent condition, 
the new bath spectral-function is calculated. With the 
updated pseudo-particle Green's functions and the new 
bath spectral-function one determines next approxima- 
tion for the T-matrix, pseudo-particle self-energies and 
local Green's function. The iteration is continued, until 
the convergence is found to the desired level. 
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C. Analytic continuation of the self-energy 



approximated the sum in Eq. (|43|) by 



The QMC simulation produces Green's function G(t) 
of imaginary time r = it or equivalently Green's func- 
tion and the self-energy defined in Matsubara frequency 
points. However, real-frequency self-energy is needed to 
obtain transport quantities. The analytic continuation 
of QMC data is required, which is an ill-posed problem 
and altogether hopeless if the precision of data is not ex- 
tremely good and if the statistical errors are not taken 
into account properly. As is well known, Pade method is 
not very useful for analytic continuation of noisy QMC 
data. The maximum entropy method (MEM)2S tries to 
overcome this problem by adding an entropy term to the 
functional to be minimized. This is one of the best meth- 
ods present available and usually produces real-frequency 
Green's function of relatively high quality provided the 
data are carefully analyzed. We refer the reader to the 
original literature for the details^. 

However, the quasiparticle peak for realistic density of 
states can have quite reach structure since at low tem- 
perature it tries to reproduce the LDA bands around 
the Fermi-level, i.e., the spectral function approaches the 
LDA density of states contracted for the quasiparticle 
renormalization amplitude Z ', A(u>) = p(w/Z + fj,o). The 
maximum entropy method has a tendency to smear out 
this reach structure because of the entropy term. At low 
temperature, this can lead to overshooting of spectral 
function and subsequently to the non-physical self-energy 
that ruins the causality. To avoid this pathology, we 
sometimes found useful to directly decompose the singu- 
lar kernel with the Singular Value Decomposition (SVD). 
When constructing the real frequency data, we took into 
account only those singular values, which are larger than 
precision of the QMC data. 

Imaginary time Green's G(r) can be expressed by the 
spectral function as 



G(r) 



or in discretized form 



diof{-uj)e- Tul A(uj) , 



(41) 



a 



SUNCA 



m< M ,r 



m>M 



QMC 



Orr 



S SUNCA 



-X^TT ,A Sl 



SUNCA 



(44) 



where M can be determined by the precision of the QMC 
data, i.e., Y^ t V t m5G t > S M - 

We plot the sum l|44(l in Fig. where first 3, 6 or 9 
coefficients were obtained from the QMC data. The cor- 
responding smallest singular value is printed in the leg- 
end of the same figure. For comparison, we also display 
the spectral function obtained by the maximum entropy 
method and the SUNCA solution for the same parame- 
ters. The difference between the various curves gives as 
a rough estimate for the accuracy of the technique. As 
we see, the quasiparticle resonance is obtained by rea- 
sonably high accuracy, while the Hubbard band is de- 
termined with less accuracy. In the inset of Fig. [5] we 
plot the same curves in a broader window. As we see, 
the SVD does not guarantee the spectra to be positive at 
higher frequencies. This however does not prevent us to 
accurately determine most of the physical quantities. 



6 0.0519 
9 0.0087 
12 0.0011 
MEM 
SUNCA 




CO/D 



to turn 

(42) 

where UW = 1 and VW = 1 are orthogonal matrices 
and S is diagonal matrix of singular values. The inversion 
is than simply given by 

A J =Y u ™4-V Tm G T . (43) 

m,r 

The magnitude of singular values drops very fast and 
only first few terms in the upper sum can be determined 
from the QMC data. The rest of the information, that 
determines mostly higher frequency points, can be ac- 
quired from the SUNCA spectral function. We therefore 



FIG. 5: Spectral function for semicircular DOS, inverse tem- 
perature j3 = 16 and density n = 0.8. Dot-dashed, full and 
double-dot-dashed curve correspond to the sum f 1441 with M 
chosen to be 6, 9 and 12, respectively. In legend, we also print 
the lowest singular value taken into account (SW). For com- 
parison we show maximum entropy spectra (dashed curve) 
and SUNCA spectra (dotted line). The inset shows the same 
spectra in a broader window. 

Within DMFT, the real frequency self-energy can be 
obtained from the local Green's function by the inversion 
of the Hilbert transform. Although the implementation 
is very straightforward, we will briefly mention the al- 
gorithm we used. In the high-frequency regime, we can 
expand the Hilbert transform in terms of moments of the 
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DOS as 



1.5 



w(z) 



D(e)de 



z — e 



E 



n+1 



(45) 



The series can be inverted and solved for z 
1 



z(w) 



w 



+ (e) + ({a') - ( e )> 



+ -3{e 2 )(e) + 2(ey)w 2 + 



(46) 



For most of the frequency points, the expansion up to 
some higher power (~ w s ) gives already an accurate es- 
timation for the inverse function. However, when w gets 
large, we need to use one of the standard root-finding 
methods to accurately determine the solution. This is 
however much easier than general root-finding in com- 
plex plane since we always have a good starting guess 
for the solution. We start evaluating the inverse function 
at high frequency where the absolute value of G is small 
and we can use the expansion in Eq. l|4^|) . Then we use 
the fact that Green's function is a continuous function of 
a real frequency and we can follow the solution from fre- 
quency point to frequency point by improving it with few 
steps of a secant (or Newton) method. A special atten- 
tion, however, must be paid not to cross the branch-cut 
and get lost in the non-physical complex plane. There- 
fore, each secant or Newton step has to be shortened if 
necessary. The self-energy is finally expressed by the in- 
verse of Hilbert transform w^ 1 as 



(47) 



Fig.Elshows the imaginary part of self-energy obtained by 
both analytic-continuation methods. As a reference and 
comparison we also show the results obtained by SUNCA 
method, which is defined and evaluated on real frequency 
axes and hence does not require analytic continuation. 
The low-frequency part of the self-energy is again very 
reliably determined and does not differ for more than 
3%. 



IV. TRANSPORT COMPUTATION 

A. Transport theory 

The transport parameters of the system are expressed 
in terms of so called kinetic coefficients, denoted here by 
A m . The equation for the electrical resistivity, p, is given 
by 



k B T 1 
e 2 A - 



(48) 



and the thermopower, S, and the thermal conductivity, 
k, are expressed through 



_kB_M 
lelAo' 



k = kf 



A Al 



(49) 
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FIG. 6: Imaginary part of the self-energy obtained from 
the Green's function by the inverse of the Hilbert transform. 
Full-line was obtained by the Singular Value decomposition, 
dashed by the maximum entropy method and dot-dashed by 
SUNCA. Parameters used are the same as in Fig. 



Within the Kubo formalism 31 the kinetic coefficients are 
given in terms of equilibrium state current-current corre- 
lation functions of the particle and the heat currents in 
the system. Namely we have 



where 



A,, 

Z (iv) 
Zx (iv) 
Z 2 (iv) 



(3 m lim Z m (iv^> uj + iO) , 



ih 
iv(3 J 



dre WT (T T f(r)f(0)) , 



— [ dTe^(T T f(T)Q*(0)) , 
WP Jo 

— f dTe^(T T Q*(r)Q*(0)) 
wp Jo 



(50) 

(51) 
(52) 
(53) 



To evaluate these correlation functions, an expression 
for electric and heat currents, j x and Q x , are needed. 
Once those currents are evaluated, then calculation of 
the transport properties within DMFT reduces to the 
evaluation of the transport function 



W = v-E Tr KWp fe ( £ K(eK( e )}: 

VP, . 



(54) 



and the transport coefficients 

A m = N spm irh decl ) xx (e)f(e)f(-e)((3er, (55) 



The momentum integral in Eq. 154f) extends over the 
Brillouin zone, Vc is the volume of the unit cell. The 
simplest form of the velocity is (kj3\— V x |feo!) = v% and 
it requires evaluation of matrix elements of V^. How- 
ever an alternative form of the current and the transport 
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function can be derived via the Peirls substitution gener- 
ally in the non-orthogonal basis and is described in Ap- 
pendix |F| These two procedures generally give different 
answers 23 i 32 i 33 . 

Next we define the energy dependent velocity as 



«fc(e) = v k - euk- 



(56) 



The second term is due to the non-orthogonality of the 
basis or more specifically due to overlap between orbitals 
at different sites, local non-orthogonality does not con- 
tribute to the velocity. The spectral density matrix p k (e) 
is the multiorbital generalization of the regular single or- 
bital density of states and is given in terms of the retarded 
Green's function, G, of the system by the equation 



2tt 



{G fc (e)-[G fc (6)]t} 



Finally the Green's function is given by 



G k (z) = [(z + n)O k 



Hi 



£(*)]' 



(57) 



(58) 



Note here, that in accordance with the DMFT the self- 
energy matrix is assumed momentum independent. Now 
given an effective Hamiltonian for the system, an over- 
lap matrix and the self-energy, the equations above give 
a complete prescription for computing the transport pa- 
rameters. For computation of Eq. 154|) we have devel- 
oped two methods, one method generalizes the analyti- 
cal tetrahedron method (ATM)2i and the other one uses 
one-particle GF method in DMFT — , used to compute 
spectral densities in band structure calculations. First 
the total Hamiltonian, H k (e) 
ized and written in the form 



Hi + E(e) is diagonal- 



H k (e) = O k A*(e)E k (e)A%(e)O k , 



(59) 



where E k is the diagonal matrix of complex eigenvalues 
and A^ and A^ are the right and the left eigenvector 
matrices respectively. Then the Green's function can be 
written as 

G k (e) = A? (e)[(e + rfl E k (e)]- l A%(e), (60) 

with / being the identity matrix. The transport function 
can now be expressed as 



2^Vc" Re ^2{ r k,qp r k,pq D k,P D k,q 



k,pq 



2 [ S k,qptt,pq + S k,pqtk,qp\ ^k,p(D kl q)*^(61) 



where the matrices r x , s x and t x are 

rl = r%{e) = A£(eK(e)A£(e), (62) 
si = 4(e)=4[( £ K(e)[4( £ )]t 
t% = t%{e) = [A%{e)]U% (e)A^(e), 

and D k is a diagonal matrix defined by 

D k =D k (e) = [(e + n)I -E k (e)]-\ (63) 



When the computation of the transport function is car- 
ried out one is faced with computing integrals of the form 



k,pq k.qp 



^ (e + fx - E kiP ){e + fi- E k<q ) ' 



(64) 



„X 4-X 

a k.pq L k 



pq lz,qp 



^ (e + /i- E kiP )(e + [i - E 



k : q! 



The strategy that is used to compute these integrals is 
similar in spirit to the analytical tetrahedron method. 
The Brillouin zone is split up into a collection of equal 
sized tetrahedra and the integral over each tetrahedron 
is taken using linear interpolation between the four cor- 
ners of the tetrahedron. In the analytical tetrahedron 
method the numerator and the energy eigenvalues in the 
denominator are linearized independently and the result- 
ing integral is then done analytically. In our case we 
would want to follow the same rule which results in two 
linear functions in the denominator. Unfortunately we 
have not been able to evaluate that integral in the most 
general case, i.e. when none of the tetrahedron corners 
are degenerate although solutions can be found for de- 
generate cases when at least two of the four corners of 
the tetrahedron are identical. Hence we have to pursue 
further approximations which we outline below. 

The two main integrals that we need to compute are 
of the form 



fceA 

TVs = E 



F(k) 



; A (z - E k>p )(z - E k , q y 
F(k) 



(65) 



keA 



{z - E k;P )(z - E kA )* 



Here A denotes the tetrahedron and SS indicates that 
the imaginary parts of both denominators have the same 
sign and OS indicates that they have the opposite sign. 
This is ensured by the fact that the self-energy is re- 
tarded and z is real. For the diagonal case (p — q) the 
Tss integral can be computed exactly by linearizing the 
eigenvalues in the denominator, one simply needs to dif- 
ferentiate the ATM formulas by Lambin and Vigneron 34 . 
For the diagonal Tos however we note that the numera- 
tor is real and therefore we can write the integral in the 
following form 



rppp 

1 OS 



LmE 



fceA 



lk,p 



1 



E, 



k,p 



(66) 



where j k .p = ImSfe^. We note that 7^ is solely due 
to the self-energy, which is momentum independent and 
thus it is reasonable to expect that j k ,p changes little 
with momentum. Hence the term in the parenthesis will 
be approximated linearly within the tetrahedron and the 
resulting integral can be computed with the ATM. 

The off-diagonal case (p ^ q) for both T55 and T Q s is 
treated the same way so we will just look at T55. Both 
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factors in the denominator are inspected and we deter- 
mine which one has larger modulus (on average if neces- 
sary). Then we write the integral on the form 



fee A 



F(k) 



1 



(z-E k ) L J [z-E k ) 



(67) 



where L indicates the denominator with the larger mod- 
ulus and S indicates the one with the smaller modulus. 
The term in the parenthesis is now approximated linearly 
within the tetrahedron and the resulting integral can be 
computed with the ATM. 

The approach described here to compute the trans- 
port function has been tested numerically against models 
where other methods can be used to evaluate the trans- 
port function. For cubic systems with nearest neighbor 
hopping one can for instance evaluate both the density 
of states and the transport function quite efficiently us- 
ing Fast Fourier Transforms!. In general the results are 
quite accurate. 



B. Small scattering limit 

In order to make connections with previous approaches 
to the computation of transport properties it is interest- 
ing to consider the small scattering limit. So we take the 
self-energy of the form 



E(e)=£'(e)+ 7 £"( e ), 



(68) 



where S'(e) is the real part of the self-energy matrix, 
7 £"(e) is the imaginary part and 7 is a small parameter. 

It is clear that the transport function will diverge as 
1 /7 and thus we can approximate the numerator matrix 
elements to zeroth order in 7. Within this approximation 
the transport function can be written as 



k,p 



>T k , p (e)S(e + fi-E ktP ), (69) 



where E 'k, P 



are the eigenvalues of H k + S'(e) and 



denotes the corresponding band velocity. The lifetime 
Tfc iP (e) is formally given by 



1 



2vr ImE, 



k,p\ 



(70) 



here E kp are the eigenvalues of the full Hamiltonian. The 
imaginary part of these eigenvalues is due to the scatter- 
ing term and is therefore to first approximation linear in 
7. The lifetime therefore diverges as 1/7 but for a finite 
value of 7 we regard Eq. (|69|l as an approximation to the 
transport function and we will refer to this approach as 
the small scattering approximation. 

In spite of the limited validity of the small scattering 
approximation it is useful in the sense that it is computa- 
tionally much simpler to evaluate the transport function 
in the small scattering approximation than in the general 



case. Therefore it can be used in order to obtain a rough 
idea of the behavior of the transport parameters. 

The equations of the small scattering approximation 
are very similar to the formulae that have been used by 
other groups to compute the transport parameters of real 
materiala^i2L2i. In particular the assumption of con- 
stant lifetime is quite often used in practice, especially 
when the thermopower is being calculated. In this case 
we obtain 



'(e) = r<F*(e), 



(71) 



where the so called transport density, <&, is defined as 

$ XX (z) = W E«p)^( £ + M - E' k>p ). (72) 



k.p 



Numerical tests have shown that while the small scat- 
tering approximation can be quite good for broad bands 
it does not work well in narrow bands such as the dynami- 
cally generated quasiparticle bands of strongly correlated 
systems due to constant time approximation used. 

In the case of the thermopower we obtain 



S = 



OO 
OO 



*o 7T 2 k B T d _ ^ TT 



$ M (e)/(e)/(-e) 



(73) 



This of course is the classical Mott relation for the ther- 
mopower. In the literature this equation is often quoted 
with the transport density replaced by the spectral den- 
sity and much emphasis placed on the fact that in case 
the Fermi-level coincides with a Van Hove singularity the 
thermopower diverges. This conclusion is not supported 
when the correct form for the thermopower is used since 
no Van Hove singularities are present in the transport 
density. 

For free electrons the transport density is given by 



$**( e ) 



and therefore we get 



12tt 2 



2m, 



3/2 



c 3/2 



s 



k B n 2 k B T 1 

"N 2 7~ F 



-2/3 



T x 0.281- 



K ' 



(74) 



(75) 



where the density, n, is measured in electrons per cu- 
bic Bohr radius and the temperature, T, is measured 
in Kelvin. In case the effective mass of the electrons is 
enhanced the thermopower will simply increase by the 
enhancement factor. 

The enhancement of the thermopower can also be de- 
duced from the Mott equation in case the only effect of 
the real part of the self-energy is to change the effective 
mass of the bands that cross the Fermi-surface. If we 
assume that the change in effective mass is the same for 
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all the bands that participate in the transport the low- 
temperature thcrmopower becomes 



ks Tr 2 kBT d 



3Z 



(76) 



where the non-interacting transport density $>°> xx (e) is 
defined by 



V c 



x\2 
P 



(77) 



Here Z denotes the quasiparticle residue of the bands 
involved. Hence we see indeed that the low-temperature 
thcrmopower is enhanced by a factor of 1/Z compared 
to the non-interacting thermopower. 

To summarize theoretical part of the paper we provide 
the computational scheme used to obtain the transport 
properties 



LDA H -^ A 



DMFT Hld A^^ TRANSPORT, 



• run LDA program to obtain the Kohn-Sham (LDA) 
Hamiltonian 

• make tight-binding calculations (or downfolding) to 
extract kinetic part of the effective Hamiltonian 

• construct the Hubbard-like effective Hamiltonian 
adding the Coulomb repulsion to the kinetic part 
of the effective Hamiltonian 

• run DMFT-QMC(SUNCA) solver to obtain the 
self-energy for the effective Hamiltonian 

• construct LDA+DMFT Hamiltonian upfolding the 
self-energy to the LDA Hamiltonian 

• run transport program with the LDA+DMFT 
Hamiltonian. 



V. TEST SYSTEM AND DMFT RESULTS 

A. Test System 

To test obtained transport equations on a realistic 
system we decided to use doped LaTiC>3 compound. 
The Lai-zSr^TiOa series has been studied very exten- 
sively in the pagtjSiiSi 4 ^* 41 ^ 4 ^* 4 ^ and can be regarded as 
being one of the prime examples exhibiting the Mott- 
Hubbard metal-insulator transition. The end compound 
LaTiC>3 when prepared well is a Mott-Hubbard insula- 
tor although in the literature it is often characterized 
as a correlated or a poor metal. At high temperature 
this material is paramagnetic. The other end compound 
SrTi03 is an uncorrelated band insulator with a direct 
gap of 3.3 eV. Electronic structure properties of the 



Lai-zSr-jTiOs series is governed by the triple degener- 
ate cubic t 2g bands of the 3d orbitals (d 1 ionic configu- 
ration)^. In the distorted structure of LaTiC>3 the de- 
generacy of the band has been lifted and the single elec- 
tron occupies a very narrow, non-degenerate d xy band 45 . 
Studies of the magnetic susceptibility do indeed indicate 
that the electronic structure of the Pbnm phase is that of 
a narrow d xy band, which then with doping changes into 
a broad t2 g band (calculated bandwidth is W = 2.7 eV) 
with degenerate d xy , d xz and d yz orbitals in the Ibmm 
and PmSm phases. As a function of doping the material 
behaves as a canonical doped Mott insulator. The spe- 
cific heat and the susceptibility are enhanced, the Hall 
number is unrenormalized while the photoemission spec- 
tral function has a resonance with a weight that decreases 
as one approaches half filling. Very near half filling, (for 
dopings less than 8 %) the physics is fairly complicated. 
At small doping an antiferromagnetic metallic phase is 
observed 42 t 46 i 47 . 

To obtain LDA band structure of LaTiC>3 we used 
the liner muffin-tin orbitals (LMTO) method in its 
atomic sphere approximation (ASA) with the basis 
Ti(4s, 4p, 3d), 0(2s,2p) and La(6s, hp, hd) assuming for 
simplicity instead a real orthorhombic structure with a 
small distortions a cubic one with the same volume and 
the lattice constant ao — 7 AO a.u.. This approximation 
brings to a slight overestimation of the effective band- 
width and underestimation of the band gap between va- 
lence and conduction bands. In photoemission studies of 
LaTiOs 12 similar basis has been used. 

Using LDA band structure one can compute and com- 
pare with experiment the linear coefficient of specific heat 
which is simply given in terms of the density of states at 
the Fermi level by 



7 = 2.357 



mj 



olK 2 



Ptot(Ef) [states/ (eV unitcell)] 



where Z is the quasiparticle residue or the inverse of the 
mass renormalization. In LDA calculations the value of 
Z is equal to one. Doping dependence of the linear coef- 
ficient of specific heat in LDA calculations was computed 
within the rigid band model. Our results along with the 
experimental data are presented in Tabled 

In general, we see that the LDA data for 7 are lower 
than the experimental values, indicating a strong mass 
renormalization. We note also that as we get closer to 
the Mott-Hubbard transition the effective mass grows 
significantly. What is consistent with DMFT modeling 
of the Mott-Hubbard transition which shows that in- 
deed the effective mass diverges at the transition. We 
should note however that this is not a necessary signa- 
ture for the Mott-Hubbard transition: in V2O3 the pres- 
sure driven metal-insulator transition is accompanied by 
the divergence of the effective mass whereas the doping 
driven transition in the same system does not show that 
divergence 4 ^ 

The physical picture of the studied material is quite 
transparent, very near half filling (dopings less than 8 %) 



Doping 


Experiment 


LDA 


5% 


16.52 


3.23 


10% 


11.51 


3.16 


20% 


8.57 


3.00 


30% 


7.70 


2.82 


40% 


6.21 


2.67 


50% 


5.38 


2.52 


60% 


4.55 


2.38 


70% 


4.35 


2.19 


80% 


3.52 


2.10 



TABLE I: The linear coefficient of specific heat, 7, for 
Lai-zSr^TiOa measured in units of — ^7^5 ■ The experimental 
data is taken from—. LDA data for the linear coefficient of 
specific heat are computed from LaTi03 LDA DOS. 



the Fermi energy becomes very small and now is compa- 
rable with the exchange interactions and structural dis- 
tortion energies. A treatment beyond single site DMFT 
then becomes important to treat the spin degrees of free- 
dom. On the other hand for moderate and large doping, 
the Kondo energy is the dominant energy and DMFT is 
expected to be accurate. This was substantiated by a 
series of papers which compared DMFT calculations in a 
single band or multiband Hubbard models using simpli- 
fied density of states with the physical properties of real 
materials. Ref. |5(] addressed the enhancement of the 
magnetic susceptibility and the specific heat as half filling 
is approached. The optical conductivity and the suppres- 
sion of the charge degrees of freedom as the Mott insula- 
tor is approached was described in Refs. l5lll52L the obser- 
vation that the Hall coefficient is not renormalized was 
found in Refs. I53ll54l Finally the thermoelectric power 
on model level using iterative perturbation theory (IPT) 
as impurity solver was investigated by Palsson et alS. 

Given the fact that only very simple tight-binding pa- 
rameterizations were used in those works, and the fact 
that a large number of experiments were fit with the 
same value of parameters one should regard the qualita- 
tive agreement with experiment as very satisfactory. The 
photoemission spectroscopy of this compound, as well as 
in other transition metal compounds are not completely 
consistent with the bulk data, and it has been argued 
that disorder, and modeling of the specific surface envi- 
ronment is required to improve the agreement with ex- 
periment — . In this situation, it is clear that this is the 
simplest system for study, and it was in fact the first 
system studied by LDA+DMFTi 

The important questions to be addressed are the de- 
gree of quantitative accuracy of DMFT. Furthermore, 
given the simplicity of this system, and the existence 
of well controlled experiments, it is an ideal system for 
testing the effects of different approximations within the 
LDA+DMFT scheme. 
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B. The model 

As we pointed out in Sec. ^ for a correct description 
of a system with strong electron correlations one needs 
to bring the self-energy into the heavy orbitals. For this 
purpose a model which correctly describes the physics 
of interacting orbitals is needed. In this paper we con- 
sider a three-band Hubbard model which underlying non- 
interacting dispersion relation is that of the degenerate 
cubic t% g band of the transition metal 3d orbitals. For 
simplicity the Hubbard interaction term is taken to be 
SU(6) invariant i.e. there is equal interaction between 
two electrons of opposite spin in the same orbital as there 
is between two electrons in different orbitals on the same 
site. The more general case will be revisited in future 
publications. 

Value of the interaction strength in our model is chosen 
large enough to exhibit metal-insulator behavior in the 
studied compound. In units of half bandwidth, D, the 
interaction strength is taken to be equal U = bD. The 
interaction strength should be regarded as an input pa- 
rameter which value has be adjusted to the experimental 
situation. Saying this, we mean the chosen interaction 
strengths should be good enough to reproduce as many 
physical properties as only possible with maximum prox- 
imity to experiment. To investigate dependence of calcu- 
lated physical properties on the interaction strength we 
calculate all properties studied in this paper for two val- 
ues of Coulomb repulsion U = 3D and 5D. On the model 
level U — AD is the value very close to the minimum in- 
teraction to get metal-insulator transition (MIT) for in- 
teger filling n — 1 even in three-fold degenerate Hubbard 
model using DMFT as an instrument which takes care of 
the interaction in the system. Hence our choice of the in- 
teraction should guarantee exploration of two physically 
different behaviours of the system with and without the 
MIT. In literature absolute value of Coulomb interaction 
is magnitude under discussion mainly because there is 
no direct and reliable method to extract it neither ex- 
perimentally no theoretically. The uncertainty between 
different theoretical method»i2iS attempted to estimate 
value of U is quite substantial ranging the interaction 
strength from 3.2 eV to 5 eV. 

It should be noticed that although this choice of pa- 
rameters is consistent with insulating behavior of this 
system it might have limited validity in the real sys- 
tem at low doping. Since Lai-^Sr^TiOs is known to 
undergo several structural transformations upon doping 
and in particular the structure of LaTiOs is distorted 
away from the cubic perovskite structure and in fact the 
distortion lifts the degeneracy of the t2 g orbitals and the 
groundstate orbital is a narrow non-degenerate d xy or- 
bital. Hence one might expect that the Mott-Hubbard 
transition in this system would be better described in a 
one-band model (x < 0.08). At larger dopings (x > 0.08) 
it is however clear that the system is degenerate and thus 
our model can be expected to give a reasonably good de- 
scription in the larger doping range. In the present paper 
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FIG. 7: LDA DOS of LaTi0 3 (dotted line with star symbols), 
its tight-binding fit(solid line) and semicircular DOS (dot- 
dashed line). Arrows indicate Fermi level position for filling 
n — 0.8 (the first one is for the semicircular DOS, and the 
second one is for the tight-binding fit). 



we do not consider the effect of lifting the degeneracy due 
to Jahn- Teller distortion rather we explore the tree-fold 
degenerate Hubbard model in the whole region of doping 
interval including n = 1 point. 

The kinetic part of the model Hamiltonian has been ob- 
tained from tight-binding LMTO ASA calculations. The 
band structure of the compound around the Fermi level 
consists of three-fold degenerate Ti 3d ti g band, host- 
ing one-electron, which is a rather well separated from 
an empty Ti 3d e g band located above ti g band. A 
rather broad gap below ti g separates Ti 3d and com- 
pletely filled 2p oxygen band. Hence it is quite straight 
forward to make the tight-binding fit of ti g band to be 
used in in the impurity solvers. To achieve asymmetry in 
tight-binding DOS one needs take into account the next 
nearest neighbors, so called, t' term on Ti sublattice. The 
dispersion which we obtained from the fit is the follow- 
ing: e k = 
where t 



2t(cos k x + cos k y )+2t' cos(k x + k y ) + 2t± cos k z , 
= -0.02424, t' = -0.006, t± = -0.00151 in Ry 
units, tig part of LaTiC>3 DOS (dotted line) and its fit 
(solid line) are presented in Fig. [7| We also added one 
more curve in Fig. [7| corresponding to semicircular DOS 
which we will use for different kind of benchmarking of 
our approach. 

Making tight-binding fit of Ti 3d ti g bands in LaTiOs 
we effectively do downfolding of the whole Hamiltonian 
(better to say the main part of the Hamiltonian without 
Ti 3d tig bands) of the compound onto three ti g bands 
i.e. we incorporate information of the whole Hamiltonian 
into part of it. In this case H e f f and O e f / are 3x3 matri- 
ces. Hhh and Ohh are parts of the Hamiltonian and the 
overlap matrix corresponding to ti g block of the original 
Hamiltonian. Hu and On are parts of the Hamiltonian 



and the overlap matrices of the Hamiltonian without ti g 
block. 

Our results of downfolding are presented in Fig. [5] 
We plotted bands of the whole Hamiltonian which corre- 
spond to Ti 3d tig bands by solid lines (the upper line is 
double degenerate one) and bands obtained from Eq. I|21|) 
by dashed lines. We see that bandwidth of ti g bands ob- 
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FIG. 8: t2 3 bands of real band structure of LaTiOa 
lines) and downfolded bands obtained from Eq. (12 1 1 
Fermi level is at ef = 12.54 eV. 
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tained from Eqs. (|21|) is only slightly smaller than the 
original one. The downfolding reproduces original band 
structure one-to-one in approximately 1 eV energy win- 
dow around the chemical potential. If one uses only Hhh 
and Ohh matrices (not effective -ff e // and O e // but sim- 
ply parts of the total Hamiltonian corresponding to ti g 
block) one gets nearly dispersionless levels only as it is 
naturally expected. 

As the last step in this section we calculate the trans- 
port function and the thermoelectrical power for three- 
fold degenerate tight-binding model. The transport func- 
tion that we obtain for this model is shown in Fig. El We 
note that the slope of the transport function in the hole 
doping (/i < 0.0) regime is positive and therefore the 
LDA evaluation of the thermopower results in electron- 
like thermopowcr. The results for a few chosen doping 
values computed at temperature equal to 300 K are dis- 
played in table [H] where we have also displayed the ex- 
perimental results from Ref. l43l 

It is quite noteworthy that for the two lowest doping 
values in the table the LDA and the experiment are in a 
good agreement. For higher values of doping however the 
experimental values are about twice as large as the LDA 
values. The good agreement at low doping should be 
regarded as mostly accidental one since the experimental 
data for doping values less than 5% show the hole-like 
thermopower which the LDA of course will not be able 
to reproduce. 
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FIG. 9: The transport function for the LDA tight-binding 
Hamiltonian computed within the small scattering approx- 
imation with a constant lifetime for all orbitals along with 
DOS for n = 0.8. Zero energy corresponds to the zero tem- 
perature Fermi-level. 




FIG. 10: Tight-binding (dashed lines) and renormalized (solid 
line) DOSes for filling n=0.8 and f3 = 32 at U = 5D. 



Doping 


Experiment 


LDA data 


5% 


-5.2 


-5.6 


25% 


-9.3 


-7.8 


50% 


-18.3 


-9.3 


75% 


-29.4 


-18.2 


80% 


-41.2 


-22.8 



TABLE II: The thermopower, S, of Lai-^Sr^TiOa at 300K 
measured in units of fiV /K is computed using LDA band 
structure. The experimental data are taken from, 43 . 



C. Summary of DMFT results 

In the previous section we described how to obtain the 
Hubbard like Hamiltonian with kinetic part coming from 
downfolded bands and interaction part defined by renor- 
malized Coulomb repulsion. In this section we address 
the issue how to compute and what is the main effect of 
the second part of the Hamiltonian, namely interactions, 
on physical properties of the studied system. The method 
which we used to solve the Hamiltonian is the dynamical 
mean- field theory which we described in section ITT1 

Strong correlations dramatically renormalize tight- 
banding LDA DOS which we plot in Fig. by dashed 
line. As it is seen from the plot it has bandwidth 2D and 
very asymmetric shape. In the same figure it is shown by 
solid line the renormalized DOS obtained from DMFT- 
QMC calculations at temperature T = 980 K. Instead 
one QP peak of the non-interacting t 2g DOS we have 
nearly classical Hubbard DOS picture where the central 
peak of the non-renormalized DOS is redistributed be- 
tween two Hubbard subbands and strongly renormalized 
central peak. Although the total bandwidth increases 



substantially width of the central peak gets strongly 
reduced. We also notice the steepness of renormalized 
DOS at the chemical potential indicating reduction of 
the quasiparticle residue which is in LDA equal to one. 

So, the main effect expected from electron interactions 
is to reproduce the Mott transition when the system ap- 
proaches an integer filling. One can sec indications of the 
MIT in filling, n, dependencies of the chemical potential, 

and quasiparticle residue, Z. The MIT is clearly in- 
dicated by jump of /j, versus n dependence (the chemical 
potential changes while filling remains the same) plotted 
in Fig. 1111 and also by vanishing energy scale seen in Z 
versus n dependence in Fig. 1121 while approaching the 
Mott transition. 

In Fig. ^] we plot the chemical potential against filling 
around filing n — 1 for three values of Coulomb inter- 
action U = 3D, AD and 5D in units of the half band- 
width and for two shapes of DOS (semicircular and tight- 
binding). We notice here that both semicircular and re- 
alistic DOS are renormalized in such a way that they run 
in interval [— D, D] with norm equal to one. First two 
upper curves presented in Fig. II ll correspond to U = 3D. 
The upper curve obtained using tight-binding DOS and 
lower one comes from semicircular DOS. The first curve 
is nearly a straight line crossing n — 1 point while the 
line corresponding to semicircular DOS is about to make 
a jump which is clearly presented in the behaviour of 
U = AD line. The jump becomes even more pronounced 
for U — 5D and both, semicircular and tight-binding, 
DOSes. Let us notice that absolute value of the jump for 
tight-binding DOS is smaller that for semicircular DOS. 
From this figure one can easily conclude that the critical 
interaction when insulating behaviour appears in the sys- 
tem should be somewhere between U = 3D and U — AD 
closer to the second value (the final conclusion about the 
insulating behaviour one can make from energy depen- 
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FIG. 11: The chemical potential, fx, versus filling, n for semi- 
circular (sc) and tight-binding (tb) DOSes and various values 
of interaction, U, and temperature f3 = 16. 
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FIG. 13: DOS at the Fermi level, p[E F ) [states /(eVunitcell)}, 
vs filling, n, for semicircular (sc) and tight-binding (tb) DOSes 
and various values of interaction, U. All of the data was 
computed for ji — 16. 
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FIG. 12: Density dependence of the quasiparticle residue, Z, 
for semicircular (sc) and tight-binding (tb) DOSes and various 
values of interaction, U and temperature (5 = 16. 



dence of DOS on real axis). 

Similar to the /i against n dependence the MIT be- 
haviour one can observe from the doping dependence of 
the quasiparticle residue which is presented in Fig. 1121 In 
Fig. ^] we plot five curves for the same values of interac- 
tion and shapes of DOSes as in the previous graph. As 
we expected for U = 3 (both DOSes) and n = 1 we have 
a finite value of Z. Notice that again (as in the previous 
plot) tight-binding DOS shows more metallic behaviour 
(larger value of Z and more straight line than in the case 
of semicircular DOS). All other values of the interaction 
clearly show insulating behaviour of the system. 

So, now when one can see how electron correlations 
change physical behaviour of the system, we remind that 
the main input to DMFT-QMC or DMFT-SUNCA pro- 



cedure consists from the shape of DOS (semicircle or 
tight-binding) and the value of interaction, U. We will 
analyze both shapes of DOSes for two values of U men- 
tioned above. 

Using our results presented in Fig. 1121 and behaviour 
of p(Ep) presented in Fig.^l we can calculate the linear 
coefficient of specific heat, 7. As we saw above, LDA 
results differ a lot from experimental values for 7. Now 
we want to know whether we can get any improvements 
applying DMFT, which changes the quasiparticle residue, 
Z, and renormalizes DOS. 

In Fig. E| we plot the linear coefficient of specific heat 
against filling for different values of Coulomb repulsion 
and semicircular and tight-binding DOSes. We can no- 
tice that for the same repulsion strength the linear coef- 
ficient of specific heat for the semicircular DOS is larger 
than for tight-binding DOS which follows from larger pin- 
ning value in the case of semicircular DOS. Comparing 
U dependencies for the semicircular DOSes we see that 
the linear coefficient of specific heat for U — 3D is nearly 
linear function till filling n = 1 which one can explain 
by almost linear dependence of the quasiparticle residue. 
For U = AD and U = 5D doping dependence of the linear 
coefficient of specific heat reproduces the experimental 
behaviour and the only question left is how close theo- 
retical and experimental results are. From the plot we 
see that in general results the semicircular DOSes are far 
from the experiment while for the tight-binding DOS the 
experimental curve just in between U — 3D and U = 5D 
lines. We can claim a rather good agreement (contrary to 
LDA situation) between DMFT and experimental curves 
for whole range of dopings. The divergence of the lin- 
ear coefficient of specific heat shows the d-electron ef- 
fective mass at the Fermi level is strongly enhanced on 
approaching the MIT—. In the case of U = 5D a small 
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"overshooting" of the linear coefficient of specific heat 
for large doping can be explained by 10-15 % inaccuracy 
in the procedure of the quasiparticle extraction (we de- 
fine it from the self-energy on Matsubara axis), plus one 
can slightly tune the interaction strength, which probably 
should be smaller something like 4.5 — 4.8D. In general 
agreement between DMFT-QMC results and experimen- 
tal one is quite good. 
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FIG. 14: The linear coefficient of specific heat, 7 [rnJ/molK 2 ], 
vs. the density for different interaction strength and DOSes 
at temperature /3 = 16. 

So we can summarize the linear coefficient of specific 
heat results saying that changes in the spectral weight, 
Z, are the main source of the improvement of our re- 
sults for the linear coefficient of the specific heat. Those 
changes are most remarkable for U — 5D where Z tends 
to zero while density approaches the integer filling n = 1. 
Diverging behaviour of linear coefficient of specific heat 
for small doping in the real material can be explained by 
one of the structural transitions happening in LaTiC>3, 
at doping less than 5% the three-fold degeneracy is lifted 
and we effectively have only one-band model for which 
U = 3D could be large enough to get the MIT transition 
at integer filling. 

Besides the degeneracy and the structural transition, 
temperature belongs to major players on the field of the 
MIT. The energy scale which separates low and high 
temperatures is given by D' = ZD where D is uncorre- 
cted bandwidth of the system and Z is the quasiparticle 
residue. Since the metal-insulator transition is accom- 
panied by the vanishing of the quasiparticle residue we 
see that in order to access the low-temperature phase we 
need to go to lower and lower temperatures as we get 
closer to integer filling. Hence we can not expect to be 
able to see clear evidence of the divergence of 7 for the 
range of temperatures we are working in QMC approach- 
ing integer filling. 



D. Comparison QMC and SUNCA 

In this section we analyze and compare two impurity 
solvers i.e. QMC and SUNCA which have been described 
in section IIIII Below we present results of those two 
methods and choose one of them to compute transport. 
Similar to the case of LDA calculation, in calculation the 
impurity problem we also need to make a compromise 
between speed and accuracy. It is well known that QMC 
impurity solver is very expensive but exact (the only ap- 
proximation used in QMC is the Trotter break up) while 
SUNCA is cheap method but it is based on more approx- 
imations. QMC works in imaginary time and Matsubara 
frequency domain while SUNCA works on real frequency 
axis. To compute transport properties one needs the self- 
energy on real axis. In the case of QMC it is necessary 
to make the analytical continuation using maximum en- 
tropy or singular decomposition method to get the self- 
energy on real axis as it was described in section 1111 CI 
This is the weakest point in the DMFT-QMC procedure. 
DMFT-SUNCA working on real axis has the self-energy 
immediately after the self-consistency is reached. 

As we noticed in the previous subsection the main task 
of the interaction (read the impurity solver) is to pro- 
duce the MIT at integer fillings. And one of the criteria 
of insulating behaviour in the system is vanishing quasi- 
particle weight. In Fig. 1151 we compare the quasiparti- 
cle residue, Z obtained from DMFT-QMC and DMFT- 
SUNCA methods as function of doping for U — 5D and 
realistic DOS. We see that both methods are in a good 
agreement with each other. We also provide Z versus 
n curve calculated using iterative perturbative theoryS 
(IPT) to see that all three impurity solvers produce at 
least qualitatively the same trends. 

Now we can go further and compare electron GF on 
Matsubara axis. Imaginary axis is a natural space of 
work for QMC and to compare results with SUNCA we 
used Lehmann representation connecting spectral func- 
tion on real axis with GF on imaginary axis. The repre- 
sentation is analytical and exact, hence, the comparison 
can be made without any assumptions and approxima- 
tions or uncertainties which could arise in the case of the 
analytical continuation. 

In Fig. 1161 we plot real and imaginary parts of GF on 
Matsubara axis for QMC (symbols) and SUNCA (lines) 
at different dopings and temperature (3 = 2 where one 
would expect very good agreement for the methods (as 
higher temperature as better and faster both methods 
work). And indeed it is the case for all the curves. In 
Fig-El wc plot GF and imaginary parts of the self-energy 
for lower temperature, (3 = 16 (temperature which is 
mostly used in our calculations) and for 10 and 20 % of 
dopings (they are our usual dopings used in calculations). 

As we can conclude from presented curves we have 
quite good agreement between the two methods and 
hence we can use SUNCA in our transport calculations 
where behaviour of the self-energy on real axis around 
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FIG. 15: Filling dependence of the quasiparticle residue, Z, 
and the linear coefficient of specific heat, 7, obtained from two 
impurity solvers: QMC (solid line with stars) and SUNCA 
(dashed line with circle symbols) for U = 5D, temperature 
f3 — 16. Experimental points are given by cross symbols 
and dot-dashed line is used as a guide for eye. Tight-binding 
density of states was used in the self-consistency loop of the 
DMFT procedure. For comparison we also provide Z vs n 
curve obtained with IPT method for the same parameters as 
ones used in QMC and SUNCA calculations. 



the Fermi level is very crucial for the transport prop- 
erties which are extremely sensitive to the shape, slope 
and value of the self-energy at the Fermi level. Transport 
properties become more and more sensitive to all the de- 
tails of the transport function at the Fermi energy with 
lowering temperature. Taking into account all the com- 
parisons made and calculations done we conclude that 
SUNCA is fast and accurate enough method in compar- 
ison with QMC to compute the transport properties of 
the compound. 



VI. RESULTS OF TRANSPORT 
CALCULATIONS 

A. Spectral and transport functions in the real 
system 

Before to do transport computations it is worth to 
study spectral and transport functions dependencies on 
doping and temperature. As we discussed in the previous 
section IV Dl we will use SUNCA as main method to com- 
pute transport properties (one can avoid the analytical 
continuation procedure in this case). But, at any rate, 





-0.2 
-0.4 


-0.5 
-1 

-1.5 
-2 



\t . Re G(ico) QMC 
1* . Im G(ico) QMC 
1 _ Re G(iffl) SUNCA 
[ - Im G(ico) SUNCA 


1 >'^*** :S r^oT = p^6 " " 

• 


1 1 

n=0.8 < ,j,ir— 


■ 1 1 

n=0.9 [3=16 | 


*%✓ . ImX(ifi)) QMC 

ImE(iffl) SUNCA - 


\r ; 

- v - 

1 





-0.2 
-0.4 


-0.5 
-1 

-1.5 







10 







10 



20 



r,2 



FIG. 16: Comparison of energy dependencies of imaginary 
and real parts of GF on Matsubara axis for different dop- 
ings computed using two impurity solvers: QMC (circles) and 
SUNCA (solid line) used in DMFT self-consistency procedure 
with semicircular DOS for U = 5D and temperature j3 = 2. 
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FIG. 17: In the two upper panels we compare energy depen- 
dencies of imaginary and real parts of GF on Matsubara axis 
for dopings n = 0.8 and 0.9 computed using QMC (circles) 
and SUNCA (solid line). In the low panels we plot imagi- 
nary parts of the self-energies for the same parameters as in 
the upper panels. We used semicircular DOS in DMFT self- 
consistency procedure and U = 5D at temperature /3 = 16. 



we also did calculations with QMC impurity solver and 
compare results obtained from the two impurity solvers 
and describe differences between them when they are the 
most noticeable. 

In Fig. ^3] we plotted the density of states per spin (the 
lower Hubbard band and quasiparticle peak are shown 
in the main panel and the inset shows the whole energy 
range) at filling n — 0.8 for various values of tempera- 
ture. Here temperature is measured in units of the half- 
bandwidth, D = 1.35 eV and thus the actual temperature 
range is quite large with the smallest temperature, corre- 
sponding to T = 0.05D, being around 780 K. The highest 
temperature plotted is equal to one but it is still not large 
enough to make incoherent motion in the system to be 
dominant. As we can see temperature changes are quite 
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substantial (the lower Hubbard band nearly disappeared 
and quasiparticle peak is shifted towards the upper Hub- 
bard band, indicating tendency to join the upper Hub- 
bard band and form incoherent broad bump) but they 
are still not enough to reach the incoherent motion state 
(the upper Hubbard band is changed but still it is very 
good separated from the QP peak - lower Hubbard band 
creation) . This situation is the expected one as we know 
the QP picture disappears for temperature higher then 
Coulomb repulsion, U, which is 5D in our case. Hence, 
for T > 5D one will see only incoherent motion in the sys- 
tem. Let us notice here the difference between SUNCA 
and QMC where in the last method the spectral density 
is just a single hump corresponding to purely incoher- 
ent carrier dynamics observed already for temperature 
(3 = 1, If we start from incoherent picture and low- 
ering temperature, then the incoherent hump splits up 
and the Hubbard bands start to form. For even lower 
temperature the lower Hubbard band moves completely 
below the Fermi-surface and the coherent quasiparticle 
peak appears at the Fermi-level. The lower Hubbard 
band starts to form at [3 = 4, QP peak is formed for 
(3 > 10. For temperature lower then (3 = 16 weight of 
the QP peak nearly does not change. We observe simi- 
lar behaviour of DOS in SUNCA where the shape of the 
QP and the lower Hubbard band change only slightly for 
temperatures lower then T = Q.1D. Described discrepan- 
cies on real axis between the two methods are entirely in 
the domain of the analytical continuation which is maxi- 
mum entropy method which reliably reproduces only low- 
energy part. One more interesting thing we should notice 
in Fig.UHlwhich is temperature dependence of DOS value 
at the Fermi level. When this value reaches the one of 
non-interacting DOS we say that the pinning condition is 
obeyed. The temperature when the pinning condition is 
reached is called the pinning temperature and it strongly 
depends on doping. For filling n = 0.8 as we can conclude 
□ 
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FIG. 18: Temperature dependence of DMFT density of states 
for n — 0.8 and U = 5D. Larger frequency interval is plotted 
in the insert. Energy is in units of half bandwidth, D. 



In Fig.^5]we plotted the density of states per spin for 
T = 0.05 and different values of doping. Choice of tem- 
perature was dictated by consideration that it should be 
lower than the pinning temperature for the largest filling 
presented. With doping growing the quasiparticle peak 
broadens and its spectral weight increases a lot while the 
weight of the lower Hubbard band changes a little (doping 
changes are 10-20 %). All the weight, which the QP peak 
gained, came from the upper Hubbard band (see insert in 
Fig. ^3 where larger energy interval is presented) . With 
doping increasing the system becomes less and less cor- 
related and in the limit of 100% doping Hubbard bands 
vanish and the quasiparticle peak transforms into free 
and empty tight-binding band. With doping decreasing 
the QP peak vanishes and the system becomes insulating 
for the repulsion U = 5D. As we mentioned it before, for 
three-band degenerate Hubbard model, U = 3-D, which 
is often used in description of LaTiOa, is not enough to 
reproduce the insulating behaviour, which appears only 
for U ps 4.5 - 5D. 
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FIG. 19: Doping dependence of DMFT density of states for 
T — 0.05 and U = 5D. Larger frequency interval is plotted 
in the insert. Energy is in units of half bandwidth, D. 

In Figs. 1201 and 1211 we presented dependence of imag- 
inary part (main panels) and real part (inserts) of the 
self-energy on temperature and doping for the same tem- 
peratures as in Figs. 1201 and the same dopings as in 
Fig. |^ In Fig. |201 we see nice quadratic behavior 
of the self-energy for low temperatures with the mini- 
mum at around the chemical potential (zero in our case) 
which is then rises and shifts with the temperature to 
the right-hand side. Real part of the self-energy reflects 
the quasiparticle residue, Z, and with lowering the tem- 
perature the QP residue increases and approaches the 
pinning value. The doping dependence of the imaginary 
part of the self-energy shows that the self-energy at the 
chemical potential decreases with increasing of doping. 
This is exactly what one should expect for a system close 
to free-electron state where more quadratic and smaller 
imaginary part of the self-energy is anticipated. The real 
part of the self-energy shows the same tendency with in- 
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creasing doping as in the case of the temperature depen- 
dence: the curve which crosses the Fermi level becomes 
more flat. At zero doping it should have zero derivative 
at the chemical potential signaling about Z = 1. The 
self-energy is extremely important characteristic of the 
system as it the only quantity which enters into trans- 
port calculations. Using the self-energy one computes 
the transport functions which is main ingredient of all 
transport equations. 
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FIG. 20: Temperature dependence of imaginary part of the 
self-energy for n = 0.8. In the inset real part of the self-energy 
is shown for the same temperatures. Energy is in units of half 
bandwidth, D. 
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FIG. 21: Doping dependence of imaginary part of the self- 
energy for T = 0.1. In the inset real part of the self-energy 
is shown for the same dopings. Energy is in units of half 
bandwidth, D. 



band and lower one plus the QP peak. But the most im- 
portant contribution to transport properties at low tem- 
peratures comes from energy region around the Fermi 
level. As it can be seen from Eqs. (|55|l the transport 
coefficients are entirely defined by the transport function 
integral in an energy window which depends on temper- 
ature. These equations allow at least qualitatively to 
define at least sign of the thermopower for small temper- 
atures. If the slope of transport function is uprising then 
the thermopower should be negative and for the other 
slope it should be positive. For large energy window the 
sign of the thermopower will strongly depend on shape 
and position to the chemical potential of the transport 
function. 
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FIG. 22: Temperature dependence of the transport function 
for n = 0.8. In the inset larger frequency interval is used. 
Energy is in units of half bandwidth, D. 
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FIG. 23: Doping dependence of the transport function for 
temperature T = 0.1. In the inset larger frequency interval is 
shown. Energy is in units of half bandwidth, D. 



In Figs. 1221 and l2~51 we plot temperature and density 
dependencies of the transport function for the same set 
of parameters as we used for Figs. |201 and [21 corre- 
spondingly. One can reveal similar features as in density 
of states: in the transport function behaviour one clearly 
identifies contributions coming from the upper Hubbard 



B. Transport parameters 

In Fig. |3| we plotted the transport parameters of the 
studied system for different densities against tempera- 
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ture. The transport parameters under consideration are 
the following: p denotes the electrical resistivity, k is 
the thermal conductivity, S is the thermopower and L 
is the Lorentz ratio. The resistivity behaviour, as it was 
found experimentally — and theoretically, is a quadratic 
function at relatively low-temperature interval becom- 
ing linear at higher temperatures (see Fig. I25fl . The 
quadratic temperature dependence of the electrical resis- 
tivity is reminiscent of the strong electron-electron scat- 
tering which predominate the elect ron-phonon scattering 
process. The thermal conductivity behaves as T~ 2 till 
temperatures of order 10 3 — 10 4 , which are relatively large 
temperatures—. The Lorentz number tends to constant 
value around 16-17 nWfl/K 2 indicating the character of 
the low-temperature scattering as Fermi-liquid one. The 
thermopower behaviour is little bit more complicated. At 
low temperature the thermopower linearly tends to zero. 
It is very hard for us to distinguish doping dependence 
for relatively small temperatures as all changes lay be- 
tween error bars which are in our case larger than the 
difference between lower and higher thermopower curves 
presented in the figure. The reason for large errors lays 
in a very small value of imaginary part of the self-energy 
which we have to deal with lowering the temperature and 
this situation is very challenging for the used impurity 
solvers. For higher temperatures (higher than 1000 K) 
we are curtain in the thermopower behaviour as there 
is no any problems with the self-energy determination 
in this temperature range. With increasing temperature 
we observe a local maximum in the temperature interval 
5 x 10 3 — 2 x 10 4 . We associate it with increasing the 
temperature cut-off (see Eq. I|55fl) which is large enough 
to take into account right-hand side slope of the central 
part in the transport function. Or in another words the 
local maximum in the thermopower in some way mimics 
behaviour of the transport function (the hump around 
the chemical potential). 
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FIG. 24: Temperature and doping dependencies of transport 
parameters: p denotes the electrical resistivity, k is the ther- 
mal conductivity, S is the thermopower and L is the Lorentz 
ratio. 

Analyzing Figs. 1241 and 1251 as functions of doping for a 
fixed temperature we can see that all curves behave in the 
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FIG. 25: The same transport parameters the same as in Fig 
1241 but on substantially larger temperature interval. 



way one would expect. The resistivity is growing with de- 
creasing doping as the system approaches the MIT while 
the thermal conductivity and the Lorentz number are de- 
creasing. The thermopower has a bit more complicated 
behaviour which depends on fixed temperature (which 
temperature slice we take). But generally it is growing 
with vanishing doping and for lower than 10 % doping one 
could get even positive thermopower which first becomes 
positive at temperatures around 5000 K (see n = 0.9 
thermopower curve) and then positiveness will propagate 
to smaller temperatures. Comparing our results with ex- 
perimental situation in doped LaTiOa^ 3 *^ we notice that 
majority of experiments are done for temperature less 
than 300 K which is a rather hard task to deal with for 
the reason we pointed out above. The biggest discrep- 
ancy we found for resistivity at small temperatures where 
our resistivity 3-4 times lower than the experimental one. 
While the thermopower behaviour (which we also treat as 
electronic one) is accurate within 30 % in absolute value. 
One would expect that the thermopower could become 
positive with decreasing doping in the way it is experi- 
mentally observed. We also could obtain it once we do 
a much more delicate and hard job taking into account 
the structural transition happening at doping x < 0.05 
as effectively we should have one-band model instead of 
three-fold degenerate. But this is beyond of the scope 
of the present work. Close to the MIT we have strongly 
asymmetric DOS and transport functions which in the 
case of integer filling n = 1 will produce positive sign 
of the thermopower. The reason for this is the position 
of negative slope (right-hand side) of the lower Hubbard 
band with is closer to the Fermi energy than the upper 
one and hence has dominant contribution into the trans- 
port properties of the system. 

So, it would be fair to say that our calculations we can 
catch at least semi-qualitative behaviour of the transport 
parameters. The electrical resistivity would require an 
additional treatment to get quantitatively a good agree- 
ment while the thermopower calculations deserve quanti- 
tative comparison with experiment and can be accurate 
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enough providing 20-30 % agreement with experiment. 



VII. CONCLUSION 

In the paper we proposed and implemented a new 
method for calculation of thermoelectrical properties in 
real materials. Dynamical mean-field theory was used 
to take into account strong electron interactions and 
thereby bring the self-energy into first-principal calcu- 
lations. Taking a rather generic for many strongly cor- 
related materials density of states, we obtained temper- 
ature and doping dependencies for such thermoelectric 
properties as electrical resistivity, the thermal conduc- 
tivity, the thermopower and the Lorentz ratio. 

We believe that the new method will be a powerful 
tool for the analysis of existing experimental data and 
guiding us to a proper physical understanding of thermo- 
electrical phenomena. This is especially important not 
only for correlated materials such as Mott-Hubbard in- 
sulators and high-temperature superconductors but also 
for simple materials like the noble metals which display 
thermoelectric behavior that still lacks a proper descrip- 
tion. In addition we hope this new method will aid in 
the search for new materials with better thermoelectri- 
cal performance by allowing for ab initio predictions of 
thermoelectric properties. 
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APPENDIX A: MANY BODY THEORY IN A 
NON-ORTHOGONAL BASIS 

Our starting point here is a representation of the ki- 
netic term of the Hamiltonian in an orthogonal basis, 
{\i}} and we assume that this basis is related to the non- 
orthogonal basis, {|a)} by the transformation matrix 

\i) = J2 \<*)Sai and (i\ = J2(a\S* ai = ]T S+(a\. 

a a a 

(Al) 

The Hamiltonian is now given by 

H = ^(iimctcj (A2) 

ij 

= SUa\H\P)S 0J cjc 3 = ^H aP c+cp. 

ija(3 a(3 

The last term in the equation above is a requirement that 
we place on the creation and destruction operators in the 



non-orthogonal basis and thus we find that 

c Z = '52 d i S L and c a =Y^S aj c 3 . (A3) 

i 3 

The non-orthogonality of the basis is encoded in the over- 
lap matrix, O a (j = (a\/3) and this matrix can be related 
to the transformation matrix, S in the following manner 

% = m = E S tMP)SfJj = £ S+OapSpj (A4) 

a/3 a/3 

Therefore we see that the overlap matrix is given by 

= (SS+)- 1 . (A5) 

We should note here that the creation operator c+ does 
not create a particle in the state \a) when acting on the 
vacuum, since as we see 

410) = £c+S£|0> = (A6) 

= £|/3)%s+ =£|/?)o^. 

03 

It is however worth noting that this state has unit overlap 
with the state \a) and zero overlap with all of the other 
basis-states! The commutation relationship of these op- 
erators are the same as for regular Fermi operators except 
that we get 

{c$M = £ S flj {ct,c j }S+ = S fii S+ = O-l (A7) 



Let us finally obtain the expression for the Green's func- 
tion in the non-orthogonal basis 



G aS3 (r) = -( Cq (t)c+(0)). 



(A8) 



The easiest way to calculate this Green's function is by 
looking at the Lagrangian for the system in the orthog- 
onal basis and then simply transform it into the non- 
orthogonal one. We have (summation over repeated in- 
dices implied) 



d 

= C+O a( 3 — Cf) - C+H a/3 Cf3. 



(A9) 



The free Matsubara Green's function can now be ob- 
tained by Fourier transforming the operators in the La- 
grangian and then the inverse of the Green's function, 
G^p(iu)), is simply the term multiplying c+Cfj. Thus we 
obtain 



G%(u) = [iuO - H]-l . 



(AlO) 
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The renormalized Green's function one gets as in the or- 
thogonal case by adding the self-energy to the Hamilto- 
nian and thus 



G(iu) = [iwO - H - E] . 



(All) 



We should remark here that these Green's functions do 
not share the same properties as their cousins in the or- 
thogonal bases do and in particular the total density is 
not given by trace of G(t = 0~). To see that we go back 
to the orthogonal basis where we know how things work 
and write the density operator as 



9 



\ ^ 



i iaf3 a/3 

(A12) 

Thus the total density of electrons in the system is 
n to t = (p) = £ O a p{c+cp) = Y G 0a(T = 0~)O a p. 

(A13) 

We should note in particular that this means that there 
seems to be no good way of assigning a density to a par- 
ticular orbital in the non-orthogonal case. 



APPENDIX B: LDA HAMILTONIAN IN 
NON-ORTHOGONAL BASE 

In LDA one has to solve the known Kohn-Sham equa- 
tion 



(-V 2 + y)* k j = ekrftfitf. 



(Bl) 



The eigenfunctions are expanded in a basis set for 
example the LMTO basis Xk( r ) which is not necessary 
orthogonal as 

a 

Substituting (j_B2|) in QBIJI we obtain 
with the normalization condition 

a/3 

The density of states with a particular character a, (3 is 
defined as 



^(e)=5>£*Of^( e -e kj . 

kj 



(B3) 



The DOS of the tight-binding fit one can get from 
Eq. (|B3(1 restricting index j to the three bands near 
the Fermi level. Fig. [7\ displays the partial t 2g DOS 
of LaTi0 3 (Eq. (|B3|) with a of t 2g character). The total 



DOS one can get from Eq. 

a, (3 



9(e) 



3} summing over indexes 
<-e ki ). (B4) 



Notice that the eigenvectors are the matrix ele- 
ments (Sgj — Ag-) of the Cholesky decomposition of 
O -1 defined in Appendix lAl and that the partial DOS 
is normalized to be of one integrated over the whole fre- 
quency range. Held et al. 13 proposed to use the partial 
t2 g DOS, the low-frequency range with a recalling so as 
to normalize it to one 



n t 2g = I 

a£t 2g 



'(e)f(e) 



(B5) 



APPENDIX C: SPLINES AND FOURIER 
TRANSFORMATIONS 

1. Direct Fourier Transformation 

In QMC program we do the direct Fourier transforma- 
tion exactly, i.e. first we obtain coefficients of the cubic 
spline exploiting physical properties of GF transformed 
and then make analytical Fourier integration knowing the 
form of the splined curve. The cubic spline interpolation 
formula reads 

G(t) = a, i +6 i (r-r i )+c i (T-T l ) 2 +d l (T-T l ) 3 , t € [n,T i+1 ], 

where coefficients cii,bi,Ci,di are equal to values of the 
function, its first, second and third derivatives at knot i 
i.e. a t = G(n), bi = G'{n), Cl - G"(n), d, = G'"{n). 

Or in terms of GF values, G; = Gfa), and its second 
derivative, = G"(ri), only 



h = 



Gi, 

Gi+i — Gi 
h 

Mj_ 

2 ' 
M l+l 



(Cl) 



2Mi + M, 



i+i 



(i 



h, 



M, 



6h 

From equations above we see that one needs to know 
the second derivatives, Mj, using tabulated values of GF, 
Gi, in order to get the cubic spline interpolation. To 
obtain Mj coefficients we use conditions of smoothness 
of the first derivative and continuity of the second one. 
As the result we have L + l equations for L + 3 unknowns 



2 A 
Mi 2 Ai 

9-2 ■ 



2 9n-l 
Pn 2 





"Mo" 




d 




Mi 




di 




.M n _ 




. d„ _ 



(C2) 
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where L is a number of time slices. In addition to L + 1 
Mo, ■■■M n , n — 0, ...L, unknowns d and rf„ also should 
be provided. Last two unknowns entirely depend on the 
boundary conditions which we have to specify in order to 
have a unique solution of Eq. (|C2|) . If one knows the first 
derivatives at the end-points then do and d n are defined 
through 



Ao - 1, d = | [^-^ - GTo) > 

MO - 1, dn-j; IG„ 



HI" 



-Gi-i 



A, 



for i G 



and d t = Ti y R ?r 

[1, n — 1]. More detailed derivations of the above formulae 
one can find in Ref. 60. 

We can reduce number of unknowns just putting M 
and M n to zero (it is so called natural spline boundary 
conditions). In this case 

A = 0, do = 0, /!„ = 0, d n — 0, 

and we have the number of unknowns matching the num- 
ber of equations, L + 1. 

This boundary condition is good enough to compute 
FT of GF in the system at or close to half filling since 
the second derivative of the Green's function is small in 
absolute value in this regime. And using the natural 
spline boundary condition we do not impose a noticeable 
error. However, away from half filling when the asym- 
metry of the system grows, along with amplitude, of one 
out of the two second derivatives, usage of the natural 
spline eventually leads to pathological behavior of the 
self-energy. The signature of this pathology is in the 
"overshooting" effect — when the self-energy at some fi- 
nite Matsubara frequency i.e. the imaginary part of the 
self-energy, becomes positive in some frequency region on 
the positive Matsubara half-axis while it should be al- 
ways negative. This, of course, amounts to having nega- 
tive spectral weight for the self-energy which is something 
that does not occur for fermionic response functions. The 
"overshooting" can get especially severe in the limiting 
cases of small temperatures, small particle densities or 
large interaction strength. 

So, to avoid the problem with the self-energy and, 
hence, with the whole procedure of the self-consistency in 
DMFT-QMC program we need to use the proper bound- 
ary conditions. And in this case we have two possibilities 
to get unique solution for the system of Eq. (|G2|) exploit- 
ing physical properties of studied GF: a) we can provide 
the first derivatives at both ends separately (in the next 
section we show how to calculate those derivatives) or b) 
we can provide the sum of the first and the sum of the 
second derivatives at the end-points, so called the first 
and the second moments of GF. 

With the second choice of the boundary conditions (b) 
the system of equations become three-diagonal one with 



two off-diagonal elements in the opposite corners of the 
matrix (— M n _i and — \M ) 



AMq + Mi 

±M +2Mi +±M 2 

\M X +2M 2 + ±M 3 
±M 2 +2M 3 



-M„_ 



= d a 
= d x 
= d 2 
= d 3 



-|M 



±M„_ 3 +2M„_ 2 iM„_i =d„- 2 

+ |M„_ 2 +2M„_! = d n _ 1 



(C3) 



where d = f + G "~°"- 1 - M«) + 2M< 2 ), 

dn-! = f ( ^-r 26 '" 1 ) - \B, G' + G' n = MM, 

AI +M n = MW. 

Now having the second derivatives Mi and, hence co- 
efficients ai,bi,Ci, di we can take Fourier integral analyt- 
ically 

G m K) = (C4) 
dr[a + b{r - r m ) + c{t - r m f - d(r - T m f)]^ lT ^ = 

T m -1 

e i-r mWn (_ 6d + 2icuj n + buj 2 - iau)3) 

- T (e iTm - lW "(-6d + 2icijj n - 6iArduj n + bu? n - 2cAtujI + 
3(AT) 2 dujl ~ iwJ n + ibArul - ^(Ar) 2 ^ + ^Ar) 3 ^)). 

Sum G m {ui n ) over m 

L 

G(uj n ) = G m (u n ), 

m— 1 

will give us the Fourier integral in frequency space. 

2. Inverse Fourier transformation 

As it is well known Green's function G(ui) falls off as 
when oj — > oo. In the program we deal with finite 
number of frequency points and cutting off 1/oj tail one 
would make a rather crude approximation as the discon- 
tinuity of GF G(t) (imaginary time domain!) has been 
removed. In such situation, the high-frequency tail has 
to be extracted from GF G(lo) and Fourier transformed 
analytically using the following Fourier relation 



1 



-[0(r)+Cn(c)]. 



(C5) 



where n(e) = l/[exp{/3e} — C] and ^ = ±1 depending on 
whether to n is bosonic or fermionic. 

The inverse Fourier transformation for GF without the 
tail is made by straightforward summation over Matsub- 
ara frequencies. Once it has been done we add the infor- 
mation about the tail using Eq. (|G5p . 



20 



APPENDIX D: MOMENTS 

Moments, M^ k \ are nothing else as the expansion of 
GF in frequency domain 



a \ V* M(k) 



(Dl) 



fc=0 



Another definition of fc-degree moment is the following 

+00 

= / dujuj k p(uj), (D2) 



where p{w) is density of states (DOS). 

Moments M^ k ' can be bind to sum of GFs and sum of 
its derivatives in imaginary-time space as 

{-l) k+1 {G {k) {0 + ) + G (k) (f3-)) = M<*\ (D3) 

where k — 0, . . . N . 

To show this one needs to take Fourier integral in parts 







G(iuj n ) = / e lu " T G{r)dT 



(D4) 



_ A (-l) fc+1 (G (fc) (0 + ) + g (fc) (/?')) 

+ { - 1)N+1 L— dN+lG{T) dr 


So, to solve the system of Eq. I|C3|I we need to adhere 
to the proper boundary conditions which are expressed 
through the various moments of the Green's function. 
What we need finally it to provide the first three moments 
M( ',MW,¥( 2 ). The first moment for Green's function 
is equal to one, the second moment proportional to the 
chemical potential in the system and the third one is a 
little bit more complicated and contains a density-density 
correlator. To show that we start with the single impurity 
Anderson model which reads 

HsiAM = ^£k a c\ a C ka +^2(e a + - ^ Ua> a )flf a 
ka a a'^ct 

+ Y. Vk ^c ka +c{j a ) (D5) 

hot 

+ y^y^ U aa '(n a n Q/ - ~(n a + n a >)), 



where e a = e a + \ J2 a ^a' U a >a, the first three moments 
be obtained from the following commutators 

mW = ({£<%;/!}+>, 



where CO — [0,H] denotes the commutator of operator 
O with the Hamiltonian, and {••■}+ is the anticommuta- 
tor. After some algebra one finds the following expres- 
sions for the moments 

M<°) = <{/„,£}> = 1, (I 
M« = ({[/ a ,fl],/t })=g a + U aa >(n a , -i), 

M (2) = ({[[f a ,H],H],fl}) = ({[f cn H},[H,fl}}) = 
(e 2 a + 2e a Y u ot<x< (n a > - -)+ 

a',a"^a k 

where J2 V ka = ~ (M3 ) 2 > anc ^ the moments Mq are 

k 

defined by Eq. (|D2|) with p(uj) = D(u), D{uj) is non- 
interacting DOS. 

Summing up similar terms in SU(N) approximation 
we get 

M (1) = e a + (2N - l)Un, (D7) 
M (2) = e 2 a + 2e a (2N- l)Un+ (D8) 

U 2 [{2N-\)n+(nn)]+Y. V L, 

k 

where n is filling per band and per spin, n = jjy nc " 
and double occupancy is defined as (nn) = ^ ("aV)' 

The second way to make the correct cubic spline as 
we mentioned before in section TG H is to provide the first 
derivatives at both ends of imaginary time interval (the 
boundary conditions). To find the first derivatives at 
the ends one can use the following definition of the first 
derivatives of finite-temperature GF 

~ <2V/„(r)/t(0)) = - (T[H, f a ]fl) = G' a (Q+). 

Using as the Hamiltonian H = Hsiam we can easily 
obtain the derivatives at the ends 

G' Q (0+) = e a {\ - n a ) + V ka c ka flj + 
^ U aa > (n a > - (n a ,n a )) , 



G' a {f3 ) = s a n a + ( ^ VkafgCka 
\ k 1 



(D9) 



where averages e.g. {J^VkaCkaflJ can be calculated 
from the following expression 

( ) = -T^2 A a (iu n )G a (iuj n ). 
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In obtained formulae (Eqs. I|D6|) - (|D9() ) we should 
know filling, n a , for each band and spin as well as we 
should know density-density correlator (n a n a >). The fill- 
ing we can extract from GF itself. Calculation of the 
correlator in QMC highlights one of the advantages of 
the method: the correlator is provided by the QMC it- 
self and one does not need to rely on any additional ap- 
proximations to obtain it as e.g. in the case of multiband 
IPT method^ where coherent potential approximation is 
used to get the correlator. At each time slice the density- 
density correlator is also computed from GF but in imag- 
inary time domain where it is simply a product of two 
Green's functions in (r, r') space. We should note here 
that we compute the correlator along with other parame- 
ters in the system at each iteration step and once the self- 
consistency is reached we have correctly obtained all the 
components and parameters in the system. And finally, 
with small enough imaginary time step At one can com- 
pletely avoid the "overshooting" problem. Keeping in 
mind the main limitation of QMC procedure U Ar/2 < 1. 



In the present computations we choose At = 1/4 which 
is good enough for the range of parameters we use in the 
current paper. 



APPENDIX E: SUNCA EQUATIONS 

In this Appendix, we explicitly give the SUNCA equa- 
tions for degenerate Anderson impurity model with N/2 
bands considering fluctuations between states with M— 1, 
M and M + 1 electrons on the impurity. Self-energies 
are analytically continued to real frequencies and pro- 
jected onto the physical Q = 1 subspace. We first define 
the ladder vertex functions T a , T b with rungs of pseudo- 
particles a (with M + 1 electrons) and b (with M — 1 
electrons), respectively, as shown diagrammatically in 
Fig. El These vertex functions, obey the following Bethe- 
Salpeter equations, 



r (w, fi) = 1 + (N — M) J def(e - Q)A°(e - n)G/(e)G (e + u- fi)T a (e, 0), (El) 
T b (u),to) = l + M I def(e-n)A° c (n-e)G f (e)G b (e + u;-n)T b (e,n), (E2) 



where /(e) is the Fermi function and A®(e) = ——ImG®(e) is the bare conduction electron density of states. For 
concreteness, all propagators are to be understood as the retarded ones. The auxiliary particle self-energies (Fig. |3J) 
are then given by, 



E/(w) = Ml def(e-Lj)A° c (Lj-e)G b (e)T a (LJ,e) 2 + (N-M) I def(e - tu)A Q c (e - u)G a {e)T b {u, e) 2 

- 2M(N — M) Jdef{e-u)Al{u-e)G b {e)Jde'f{e'-e)Al{e'-e)Gf{e')G a {e' + w-e), (E3) 
Sb(cj) = (N-M + l) J def(e-u)A°(e-u>)Gf(e)T a (e,u)) + (N-M + i)(N-M) J def(e - cj)A° c (c - w)G f (e) x 
de'f(e' - w)A° e (f? - uj)G f (e')G a (e' + e - u) [T b (e, e' + e-u) T b {e' , e' + e - u) - 1] , (E4) 
S a H = (A/ + 1) / def(e-u})A° c (w-e)Gf(e)T b (e,Lj) + {M + 1)M J def(e-uj)A c (u-e)G f (e)x 

de'f(e' - lo)A° c (lu - e')G f (e')G b (e' + e - w) [T a (e, e' + e-uj) T a (e', e' + e - u) - 1] . (E5) 



In order to calculate the local electron spectral function Ad from the self-consistently determined G a , G b , Gf, it is 
convenient to define modified vertex functions as 

SfCw.fi) = 1 + (N-M) J def(e-n)A° c (e-n)Re{G f (e)T a (e,n)}G a (e + to), (E6) 

Si(fiJ,Q) = (N-M) J def{e-n)A Q c (e-n)ln,{G f (e)T a {e,ti)}G a {e + Lo), (E7) 

S?(w,ft) = l + M f def(e-n)A c (n-e)Re{Gf(e)T b (e,n)}G b (e-u), (E8) 

Sl(u,Sl) = M J def(e-n)A° c {n-e)Im{G f {e)T b {e,ty}G b (e-uj). (E9) 
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The local spectral function then reads 
A d {uo) = - 



M-i ) Im J ^jiz^ff( n + "HMGb(ms«(u,n) 2 -si(u,n) 2 }+ 2Re[G b (n)]s?(u,n)s I a (u;,n)} 



AT — 1 \ f dO p~@ n 

M J Im 7 ^ 7R G/( " ~ W) {MG.fflSf (w,fi) 2 ~ #(w,n) 2 ]+ 2Re[G (fi)]Sf (u;,n)Sfta;,n)} 



+ 2M 



N - 1 
M 



def(e - fl)A° c (e - tt) Im[G b (Q)Gf(e)]Im[Gf(Q + uj)G a (e + u)]. 



(E10) 



Note that the exponential divergencies of the statistical 
factors appearing in Eq. (|E10() are compensated by the 
threshold behavior of the corresponding auxiliary parti- 
cle spectral functions A M (w) = — —Im.G ^(lj) , fi = a,b,f 
in the integrands. For the numerical treatment, these 
divergencies can be explicitly absorbed by formulating 
the self-consistency equations (Al)-(AIO) in terms of the 
functions A^(jx>) which are defined via 



AJw) = f(-uj)A^), 



(Ell) 



and, hence, have no exponential divergence. We thus 
have, e.g., exp(-/3w)A M (w) = 



k — A p and k — > k + Ah d T respectively. This replace- 
ment is performed in both the kinetic and the interaction 
terms but not in the field operators. In our case however 
the overlap matrix appearing in the action depends also 
on momentum and therefore the proper generalization of 
the currents to non-orthogonal basis will also take the 
overlap matrix into account. Thus we obtain 



Q = 



d{Q[A p ] Or +H[A P }) 
dl p 

d{Q[A h ] d T +H[A h ]) i 



u„=o> 



lA h =0- 



(F3) 
(F4) 



APPENDIX F: TRANSPORT CALCULATIONS: 
CURRENTS DERIVATION 



Performing these operations leads to the following ex- 
pressions 



Below we derive the expressions for the currents in 
a general basis. This is done by extending the gauge- 
theoretic method developed in Ref. Jj In the non- 
orthogonal basis the action for the system can be ex- 
pressed as follows 



/ dr c\ a (o kaf 3 d T +H ka pj c k p. (Fl) 

^ k 



Here dr = \ {d T — d T ) denotes the anti-symmetrized time 
derivative. The particle and heat currents can now be 
obtained by considering the invariance of the action un- 
der local phase transformation and local translations in 
time respectively. In the orthogonal case one is lead to 
the following expression for the currents 



dH\A„ 



J = 



dA t 



and Q = — 



dH[A h ] 



A„=0 



dA t 



A h =0 



(F2) 

where A p and Ah are gauge fields conjugate to the cur- 
rents and H[Ap] and H[Ah] denote the gauged Hamil- 
tonian, i.e. the Hamiltonian with the replacements k — > 



(^a/3 S l°i/3 - *k,afi B kifi 



kail 



ka/3 

where we have defined 



B 



and 



(F5) 
(F6) 

(F7) 
(F8) 



where H® a g is the tight-binding, LMTO Hamiltonian 
of the system and O k ,a/3 is the overlap matrix that cap- 
tures the non-orthogonality of the basis that we are using. 
The validity of the expressions above is not restricted to 
DMFT and they are in fact true for all density-density 
interactions such as the Hubbard interaction. This is 
because the interaction terms are gauge invariant and 
therefore they do not contribute to the expressions for 
the currents. 
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